GeoSubSampler.py

# -*- coding: utf-8 -*-
"""
/***************************************************************************
 GeoSubSampler
                                 A QGIS plugin
 Subsamples geological orientations, faults and polygons
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2025-07-28
        git sha              : $Format:%H$
        copyright            : (C) 2025 by Ranee Joshi
        email                : raneejoshi@gmail.com
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""
from qgis.PyQt.QtCore import QSettings, QTranslator, QCoreApplication, Qt
from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtWidgets import QAction
from qgis.core import QgsProject, QgsMapLayerProxyModel, Qgis
from qgis.core import QgsVectorLayer, QgsProject
from qgis.PyQt.QtWidgets import QDockWidget
from qgis.core import QgsMessageLog, Qgis


# Qt5/Qt6 Compatibility Layer
try:
    # Try Qt6 style first
    _test = Qt.DockWidgetArea.RightDockWidgetArea
    # Qt6 detected
    QT6 = True

    # Qt6 style enums are already available
    RightDockWidgetArea = Qt.DockWidgetArea.RightDockWidgetArea


except AttributeError:
    # Qt5 detected
    QT6 = False

    # Qt5 style enums
    RightDockWidgetArea = Qt.RightDockWidgetArea


# Initialize Qt resources from file resources.py
from .resources import *

# Import the code for the DockWidget
from .GeoSubSampler_dockwidget import GeoSubSamplerDockWidget
import os.path
import tempfile
import shutil
from .calcs.StructuralPolygonSubSampler import StructuralPolygonSubSampler
from .calcs.PolygonTriangulator import PolygonTriangulator
from .calcs.FaultLineMerger import FaultLineMerger
from .calcs.FaultLengths import FaultLengths
from .calcs.FaultsGraph import FaultsGraph
from .calcs.FaultStratOffset import FaultStratOffset
from .calcs.FaultClusterOrientation import FaultsOrientations
from .calcs.FirstOrderOrientation import SubsamplingEngine, save_grid_to_shapefile
from .calcs.TopferPillewizer import topfer_count, scale_iterations, scale_points_tp, scale_lines_tp, scale_polygons_tp
from .calcs.PolygonSimplification import SimplificationEngine, vector_simplify_file_two_stage
import geopandas as gpd
import os
import random
import re
import time
import math
import pandas as pd
import numpy as np


class GeoSubSampler:
    """QGIS Plugin Implementation."""

    def __init__(self, iface):
        """Constructor.

        :param iface: An interface instance that will be passed to this class
            which provides the hook by which you can manipulate the QGIS
            application at run time.
        :type iface: QgsInterface
        """
        # Save reference to the QGIS interface
        self.iface = iface

        # initialize plugin directory
        self.plugin_dir = os.path.dirname(__file__)

        # initialize locale
        locale = QSettings().value("locale/userLocale")[0:2]
        locale_path = os.path.join(
            self.plugin_dir, "i18n", "GeoSubSampler_{}.qm".format(locale)
        )

        if os.path.exists(locale_path):
            self.translator = QTranslator()
            self.translator.load(locale_path)
            QCoreApplication.installTranslator(self.translator)

        # Declare instance attributes
        self.actions = []
        self.menu = self.tr("&GeoSubSampler")
        # TODO: We are going to let the user set this up in a future iteration
        self.toolbar = self.iface.addToolBar("GeoSubSampler")
        self.toolbar.setObjectName("GeoSubSampler")

        # print "** INITIALIZING GeoSubSampler"

        self.pluginIsActive = False
        self.dockwidget = None

    # noinspection PyMethodMayBeStatic
    def tr(self, message):
        """Get the translation for a string using Qt translation API.

        We implement this ourselves since we do not inherit QObject.

        :param message: String for translation.
        :type message: str, QString

        :returns: Translated version of message.
        :rtype: QString
        """
        # noinspection PyTypeChecker,PyArgumentList,PyCallByClass
        return QCoreApplication.translate("GeoSubSampler", message)

    def add_action(
        self,
        icon_path,
        text,
        callback,
        enabled_flag=True,
        add_to_menu=True,
        add_to_toolbar=True,
        status_tip=None,
        whats_this=None,
        parent=None,
    ):
        """Add a toolbar icon to the toolbar.

        :param icon_path: Path to the icon for this action. Can be a resource
            path (e.g. ':/plugins/foo/bar.png') or a normal file system path.
        :type icon_path: str

        :param text: Text that should be shown in menu items for this action.
        :type text: str

        :param callback: Function to be called when the action is triggered.
        :type callback: function

        :param enabled_flag: A flag indicating if the action should be enabled
            by default. Defaults to True.
        :type enabled_flag: bool

        :param add_to_menu: Flag indicating whether the action should also
            be added to the menu. Defaults to True.
        :type add_to_menu: bool

        :param add_to_toolbar: Flag indicating whether the action should also
            be added to the toolbar. Defaults to True.
        :type add_to_toolbar: bool

        :param status_tip: Optional text to show in a popup when mouse pointer
            hovers over the action.
        :type status_tip: str

        :param parent: Parent widget for the new action. Defaults None.
        :type parent: QWidget

        :param whats_this: Optional text to show in the status bar when the
            mouse pointer hovers over the action.

        :returns: The action that was created. Note that the action is also
            added to self.actions list.
        :rtype: QAction
        """

        icon = QIcon(icon_path)
        action = QAction(icon, text, parent)
        action.triggered.connect(callback)
        action.setEnabled(enabled_flag)

        if status_tip is not None:
            action.setStatusTip(status_tip)

        if whats_this is not None:
            action.setWhatsThis(whats_this)

        if add_to_toolbar:
            self.toolbar.addAction(action)

        if add_to_menu:
            self.iface.addPluginToMenu(self.menu, action)

        self.actions.append(action)

        return action

    def initGui(self):
        """Create the menu entries and toolbar icons inside the QGIS GUI."""

        icon_path = ":/plugins/GeoSubSampler/icon.png"
        self.add_action(
            icon_path,
            text=self.tr("GeoSubSampler"),
            callback=self.run,
            parent=self.iface.mainWindow(),
        )

    # --------------------------------------------------------------------------

    def onClosePlugin(self):
        """Cleanup necessary items here when plugin dockwidget is closed"""

        # print "** CLOSING GeoSubSampler"

        # disconnects
        self.dockwidget.closingPlugin.disconnect(self.onClosePlugin)

        # remove this statement if dockwidget is to remain
        # for reuse if plugin is reopened
        # Commented next statement since it causes QGIS crashe
        # when closing the docked window:
        # self.dockwidget = None

        self.pluginIsActive = False

    def unload(self):
        """Removes the plugin menu item and icon from QGIS GUI."""

        # print "** UNLOAD GeoSubSampler"

        for action in self.actions:
            self.iface.removePluginMenu(self.tr("&GeoSubSampler"), action)
            self.iface.removeToolBarIcon(action)
        # remove the toolbar
        del self.toolbar

    # --------------------------------------------------------------------------

    # --------------------------------------------------------------------------
    def removeAllLayersByName(self, layerName):
        """
        Remove all layers with the given name
        """
        layers = QgsProject.instance().mapLayersByName(layerName)

        if layers:
            for layer in layers:
                QgsProject.instance().removeMapLayer(layer.id())
            # print(f"Removed {len(layers)} layer(s) named '{layerName}'")
            return len(layers)
        else:
            print(f"No layers found with name '{layerName}'")
            return 0

    @staticmethod
    def _shapefile_safe_columns(columns):
        """
        Map each column name to a shapefile/DBF-safe (<=10 char) name,
        de-duplicating collisions with a numeric suffix. Applied ourselves
        before writing, so the result is exactly known (rather than relying
        on whatever the OGR shapefile driver would silently truncate names to).
        """
        seen = {}
        mapping = {}
        for col in columns:
            if col == 'geometry':
                mapping[col] = col
                continue
            candidate = str(col)[:10]
            n = 1
            while candidate in seen and seen[candidate] != col:
                suffix = f"_{n}"
                candidate = str(col)[:10 - len(suffix)] + suffix
                n += 1
            seen[candidate] = col
            mapping[col] = candidate
        return mapping

    def _resolved_source(self, layer):
        """
        Return a .shp path for `layer`, usable with os.path.exists()/gpd.read_file().

        Every tool in this plugin assumes a shapefile-backed source. A
        GeoPackage layer's source string (e.g. "C:/foo.gpkg|layername=bar")
        fails os.path.exists() outright, so such layers were silently
        treated as missing and skipped. If the source isn't already a .shp,
        export the layer's current features to a sibling .shp next to the
        geopackage (never written back into it) and reuse that file on
        repeat runs. Field names longer than the 10-char shapefile/DBF limit
        are pre-truncated ourselves (see _shapefile_safe_columns) and the
        original->truncated mapping is cached for _resolved_field() to use.
        """
        source = layer.source()
        path_part = source.split('|')[0]

        if path_part.lower().endswith('.shp'):
            return path_part

        if not os.path.exists(path_part):
            return source

        out_dir = os.path.dirname(path_part)
        base = os.path.splitext(os.path.basename(path_part))[0]
        safe_name = re.sub(r'[^A-Za-z0-9_]+', '_', layer.name()).strip('_')
        out_path = os.path.join(out_dir, f"{base}_{safe_name}.shp")

        if not hasattr(self, '_field_name_maps'):
            self._field_name_maps = {}

        if not os.path.exists(out_path) or out_path not in self._field_name_maps:
            match = re.search(r'layername=([^|]+)', source)
            layername = match.group(1) if match else None
            try:
                gdf = (gpd.read_file(path_part, layer=layername) if layername
                       else gpd.read_file(path_part))
                field_map = self._shapefile_safe_columns(gdf.columns)
                self._field_name_maps[out_path] = field_map
                if not os.path.exists(out_path):
                    if any(k != v for k, v in field_map.items()):
                        gdf = gdf.rename(columns=field_map)
                    gdf.to_file(out_path, driver='ESRI Shapefile')
            except Exception as exc:
                print(f"Could not export layer '{layer.name()}' to shapefile: {exc}")
                return source

        return out_path

    def _resolved_field(self, layer, field_name):
        """
        Translate a field combo-box selection into the column name actually
        present in `layer`'s resolved/materialized data. Handles two
        independent name changes:

          1. QgsFieldComboBox displays a field's alias if one is set on the
             layer, not its real name — resolved via QgsFields.lookupField(),
             which is alias-aware (common on curated/compiled datasets where
             friendly aliases are configured over cryptic real column names).
          2. _resolved_source's shapefile export truncates real names longer
             than 10 characters — resolved via its cached field-name map.
        """
        if not field_name:
            return field_name

        real_name = field_name
        try:
            fields = layer.fields()
            idx = fields.lookupField(field_name)
            if idx >= 0:
                real_name = fields.at(idx).name()
        except Exception:
            pass

        out_path = self._resolved_source(layer)
        field_map = getattr(self, '_field_name_maps', {}).get(out_path)
        if field_map and real_name in field_map:
            return field_map[real_name]
        return real_name

    def setUpPointSampler(self, checkFields=True):
        layerName = self.dockwidget.mMapLayerComboBox_points.currentText()
        dip_col = self.dockwidget.mFieldComboBox_dip.currentText()
        dip_dir_col = self.dockwidget.mFieldComboBox_dip_dir.currentText()

        if layerName == "":
            return False

        if os.path.exists(self._resolved_source(self.points_layer)):

            gdf = gpd.read_file(self._resolved_source(self.points_layer))

            if (dip_col == "" or dip_dir_col == "") and checkFields:
                self.iface.messageBar().pushMessage(
                    "Please define dip & dip dir/strike fields before continuing!",
                    level=Qgis.Warning,
                    duration=15,
                )
                return False
            elif dip_col == "" or dip_dir_col == "":
                dip_col = gdf.columns[1]
                dip_dir_col = gdf.columns[2]
            else:
                # Field combo boxes may show a field's alias rather than its
                # real name, and/or the real name may have been shortened by
                # _resolved_source's shapefile export (10-char DBF limit) —
                # _resolved_field() untangles both.
                dip_col = self._resolved_field(self.points_layer, dip_col)
                dip_dir_col = self._resolved_field(self.points_layer, dip_dir_col)

            missing = [c for c in (dip_col, dip_dir_col) if c not in gdf.columns]
            if missing:
                self.iface.messageBar().pushMessage(
                    f"T&P: field(s) {missing} not found in the loaded data. "
                    f"Available columns: {list(gdf.columns)}",
                    level=Qgis.Critical, duration=20)
                return False

            # Add coordinate columns from geometry (engine needs these for grid operations)
            gdf = gdf.copy()
            gdf['EASTING'] = gdf.geometry.x
            gdf['NORTHING'] = gdf.geometry.y

            # Drop features with nulls in dip/dip-dir before they reach the engine.
            # Unfiltered nulls become NaN vector components; pandas' skipna sum then
            # silently excludes them from their grid cell's mean/Kent calculation,
            # which can make whole cells vanish from the subsampled output.
            n_before = len(gdf)
            gdf = gdf.dropna(subset=[dip_col, dip_dir_col, 'EASTING', 'NORTHING']).copy()
            n_dropped = n_before - len(gdf)
            if n_dropped:
                self.iface.messageBar().pushMessage(
                    f"{n_dropped} point(s) with null dip/dip-dir/coordinates "
                    f"dropped before subsampling.",
                    level=Qgis.Warning, duration=10)

            # If input is strike, convert to dip direction (dip_dir = strike - 90)
            if not self.dockwidget.checkBox_dip_dir.isChecked():
                gdf['_eff_dipdir'] = (gdf[dip_dir_col].astype(float) - 90) % 360
                eff_dipdir_col = '_eff_dipdir'
            else:
                eff_dipdir_col = dip_dir_col

            engine = SubsamplingEngine(
                dip=dip_col, dipdir=eff_dipdir_col,
                easting='EASTING', northing='NORTHING'
            )
            return (gdf, engine, dip_col, dip_dir_col)
        else:
            self.iface.messageBar().pushMessage(
                "Sorry, only layers saved to disk as a shapefile can be processed!",
                level=Qgis.Warning,
                duration=15,
            )
            return False

    def _add_layer_with_style(self, new_layer, source_layer, fail_msg=''):
        """Add new_layer to the project and copy style from source_layer."""
        if not new_layer.isValid():
            print(f"Failed to load layer: {fail_msg or new_layer.name()}")
            return
        try:
            from qgis.PyQt.QtXml import QDomDocument
            doc = QDomDocument()
            source_layer.exportNamedStyle(doc)
            new_layer.importNamedStyle(doc)
        except Exception:
            pass
        QgsProject.instance().addMapLayer(new_layer)

    def finalisePointSampler(self, gdf2, qgis_layer, name, param=0):
        layer_path = os.path.dirname(self._resolved_source(qgis_layer))
        if(param==0):
            param2 = ''

        else:
            param2="_"+ str(int(param))

        new_path = (
            layer_path
            + "/"
            + qgis_layer.name()
            + "_"
            + name
            + param2
            + ".shp"
        )
        # print("new_path:", new_path)
        if os.path.exists(new_path):
            random_5_digit_integer = random.randint(10000, 99999)
            new_path = (
                layer_path
                + "/"
                + qgis_layer.name()
                + "_"
                + name
                + param2
                + "_"
                + str(random_5_digit_integer)
                + ".shp"
            )
        # print("new_path_2:", new_path)
        # Write GeoPandas back to file
        gdf2.to_file(new_path, driver="ESRI Shapefile")

        # reload as layer
        upscaled_layer = QgsVectorLayer(
            new_path, qgis_layer.name() + "_" + name +  str(param2), "ogr"
        )

        self._add_layer_with_style(upscaled_layer, qgis_layer,
                                   qgis_layer.name() + "_" + name)

    def decimation(self):
        result = self.setUpPointSampler()
        if not result:
            return
        gdf, engine, dip_col, dip_dir_col = result
        n = int(self.dockwidget.lineEdit_rate.text())
        tmpdir = tempfile.mkdtemp()
        try:
            gdf2 = engine.decimation(gdf, n=n, path_out=tmpdir)
            gdf2 = gdf2.set_crs(gdf.crs, allow_override=True)
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)
        gdf2 = gdf2.drop(columns=[c for c in ['EASTING', 'NORTHING', '_eff_dipdir'] if c in gdf2.columns])
        self.finalisePointSampler(gdf2, self.points_layer, "decimate", n)

    def stochastic(self):
        result = self.setUpPointSampler(False)
        if not result:
            return
        gdf, engine, dip_col, dip_dir_col = result
        retain_percentage = float(self.dockwidget.mQgsDoubleSpinBox_percent.value())
        frac = retain_percentage / 100.0
        tmpdir = tempfile.mkdtemp()
        try:
            gdf2 = engine.stochastic(gdf, frac=frac, random_state=int(time.time()),
                                     path_out=tmpdir)
            gdf2 = gdf2.set_crs(gdf.crs, allow_override=True)
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)
        gdf2 = gdf2.drop(columns=[c for c in ['EASTING', 'NORTHING', '_eff_dipdir'] if c in gdf2.columns])
        self.finalisePointSampler(gdf2, self.points_layer, "stochastic", retain_percentage)

    def _count_gridcell_output(self, run_fn, gdf, min_x, max_x, min_y, max_y, cs):
        """Run engine grid-cell method at cell size cs and return valid point count."""
        tmpdir = tempfile.mkdtemp()
        try:
            file = run_fn(gdf, min_x, max_x, min_y, max_y, n=cs, path_out=tmpdir)
            df = pd.read_csv(os.path.join(tmpdir, file + ".csv"))
            return int((df['DIP'] != -999).sum()) if not df.empty else 0
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)

    def _find_grid_size_for_target(self, run_fn, gdf, min_x, max_x, min_y, max_y, target_n):
        """
        Find the grid cell size that produces approximately target_n output points.

        Strategy (Option 3):
          1. Analytical estimate: cell_size = sqrt(area / target_n)
          2. Proportional correction using actual vs target count
          3. Binary search if still outside 15% tolerance
        Returns (cell_size_in_crs_units, actual_count).
        """
        area = (max_x - min_x) * (max_y - min_y)
        tol = 0.15

        def count(cs):
            n = self._count_gridcell_output(run_fn, gdf, min_x, max_x, min_y, max_y, cs)
            print(f"  [target N] cell={cs:.2f} → {n} points")
            return n

        # Step 1: analytical estimate
        cs1 = math.sqrt(area / max(1, target_n))
        n1 = count(cs1)
        if n1 == target_n or abs(n1 - target_n) <= tol * target_n:
            return cs1, n1

        # Step 2: proportional correction (count ∝ 1/cell_size²)
        cs2 = cs1 * math.sqrt(max(1, n1) / max(1, target_n))
        n2 = count(cs2)
        if n2 == target_n or abs(n2 - target_n) <= tol * target_n:
            return cs2, n2

        # Step 3: binary search
        # lo (small cell) → many points, hi (large cell) → few points
        if cs1 < cs2:
            lo, lo_n, hi, hi_n = cs1, n1, cs2, n2
        else:
            lo, lo_n, hi, hi_n = cs2, n2, cs1, n1

        for _ in range(8):
            if lo_n >= target_n:
                break
            lo /= 2
            lo_n = count(lo)
        for _ in range(8):
            if hi_n <= target_n:
                break
            hi *= 2
            hi_n = count(hi)

        best_cs = cs2 if abs(n2 - target_n) < abs(n1 - target_n) else cs1
        best_n  = n2  if abs(n2 - target_n) < abs(n1 - target_n) else n1

        for _ in range(12):
            if hi - lo < 0.1:
                break
            mid = (lo + hi) / 2
            mid_n = count(mid)
            if abs(mid_n - target_n) < abs(best_n - target_n):
                best_cs, best_n = mid, mid_n
            if abs(mid_n - target_n) <= tol * target_n or mid_n == target_n:
                break
            if mid_n > target_n:
                lo = mid
            else:
                hi = mid

        return best_cs, best_n

    def gridCellAveraging(self):
        result = self.setUpPointSampler()
        if not result:
            return
        gdf, engine, dip_col, dip_dir_col = result
        self.point_layer = self.dockwidget.mMapLayerComboBox_points.currentLayer()
        is_geographic = self.point_layer.crs().isGeographic()
        bounds = gdf.total_bounds
        min_x, max_x = int(bounds[0]), int(bounds[2]) + 1
        min_y, max_y = int(bounds[1]), int(bounds[3]) + 1

        if self.dockwidget.checkBox_target_n.isChecked():
            target_n = self.dockwidget.spinBox_target_n.value()
            grid_size, _ = self._find_grid_size_for_target(
                engine.gridcell_average, gdf, min_x, max_x, min_y, max_y, target_n)
            display_size = grid_size * 110000 if is_geographic else grid_size
            self.dockwidget.lineEdit_grid_size.setText(f"{display_size:.0f}")
        else:
            grid_size = float(self.dockwidget.lineEdit_grid_size.text())
            display_size = grid_size
            if is_geographic:
                grid_size = grid_size / 110000

        tmpdir = tempfile.mkdtemp()
        try:
            file = engine.gridcell_average(gdf, min_x, max_x, min_y, max_y,
                                           n=grid_size, path_out=tmpdir)
            gdf2 = save_grid_to_shapefile(tmpdir, file)
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)
        if gdf2 is not None:
            gdf2 = gdf2.set_crs(gdf.crs, allow_override=True)
            gdf2 = gdf2.drop(columns=['EASTING', 'NORTHING'], errors='ignore').rename(
                columns={'DIP': dip_col, 'DIP_DIR': dip_dir_col})
            self.finalisePointSampler(gdf2, self.points_layer, "gridCellAveraging", display_size)

    def kent(self):
        result = self.setUpPointSampler()
        if not result:
            return
        gdf, engine, dip_col, dip_dir_col = result
        self.point_layer = self.dockwidget.mMapLayerComboBox_points.currentLayer()
        is_geographic = self.point_layer.crs().isGeographic()
        bounds = gdf.total_bounds
        min_x, max_x = int(bounds[0]), int(bounds[2]) + 1
        min_y, max_y = int(bounds[1]), int(bounds[3]) + 1

        if self.dockwidget.checkBox_target_n.isChecked():
            target_n = self.dockwidget.spinBox_target_n.value()
            grid_size, _ = self._find_grid_size_for_target(
                engine.spherical_kent, gdf, min_x, max_x, min_y, max_y, target_n)
            display_size = grid_size * 110000 if is_geographic else grid_size
            self.dockwidget.lineEdit_grid_size_kent.setText(f"{display_size:.0f}")
        else:
            grid_size = float(self.dockwidget.lineEdit_grid_size_kent.text())
            display_size = grid_size
            if is_geographic:
                grid_size = grid_size / 110000

        tmpdir = tempfile.mkdtemp()
        try:
            file = engine.spherical_kent(gdf, min_x, max_x, min_y, max_y,
                                         n=grid_size, path_out=tmpdir)
            gdf2 = save_grid_to_shapefile(tmpdir, file)
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)
        if gdf2 is not None:
            gdf2 = gdf2.set_crs(gdf.crs, allow_override=True)
            gdf2 = gdf2.drop(columns=['EASTING', 'NORTHING'], errors='ignore').rename(
                columns={'DIP': dip_col, 'DIP_DIR': dip_dir_col})
            self.finalisePointSampler(gdf2, self.points_layer, "grid_cell_kent", display_size)

    def kentOutlier(self):
        result = self.setUpPointSampler()
        if not result:
            return
        gdf, engine, dip_col, dip_dir_col = result
        self.point_layer = self.dockwidget.mMapLayerComboBox_points.currentLayer()
        is_geographic = self.point_layer.crs().isGeographic()
        bounds = gdf.total_bounds
        min_x, max_x = int(bounds[0]), int(bounds[2]) + 1
        min_y, max_y = int(bounds[1]), int(bounds[3]) + 1

        if self.dockwidget.checkBox_target_n.isChecked():
            target_n = self.dockwidget.spinBox_target_n.value()
            grid_size, _ = self._find_grid_size_for_target(
                engine.outlier_removal, gdf, min_x, max_x, min_y, max_y, target_n)
            display_size = grid_size * 110000 if is_geographic else grid_size
            self.dockwidget.lineEdit_grid_size_kent_2.setText(f"{display_size:.0f}")
        else:
            grid_size = float(self.dockwidget.lineEdit_grid_size_kent_2.text())
            display_size = grid_size
            if is_geographic:
                grid_size = grid_size / 110000

        tmpdir = tempfile.mkdtemp()
        try:
            file = engine.outlier_removal(gdf, min_x, max_x, min_y, max_y,
                                          n=grid_size, path_out=tmpdir)
            gdf2 = save_grid_to_shapefile(tmpdir, file)
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)
        if gdf2 is not None:
            gdf2 = gdf2.set_crs(gdf.crs, allow_override=True)
            gdf2 = gdf2.drop(columns=['EASTING', 'NORTHING'], errors='ignore').rename(
                columns={'DIP': dip_col, 'DIP_DIR': dip_dir_col})
            self.finalisePointSampler(gdf2, self.points_layer, "grid_cell_kentOutlier", display_size)

    def _validate_dip_fields(self, require_points_layer):
        """
        Fail fast with a clear message if dip/dip-direction fields aren't
        set, before any processing (GeoPackage materialization, baseline
        output writes, etc.) has run.

        require_points_layer=True: a points layer must also be selected
        (used by the standalone point-subsampling tools, where points are
        the whole point of the operation).
        require_points_layer=False: a points layer is optional (used by
        T&P, which can run on just faults/polygons) — dip/dip-direction are
        only required if a points layer has actually been chosen.
        """
        layer_name = self.dockwidget.mMapLayerComboBox_points.currentText()
        if not layer_name:
            if require_points_layer:
                self.iface.messageBar().pushMessage(
                    "Please select a points layer before continuing!",
                    level=Qgis.Warning, duration=15)
                return False
            return True

        dip_col = self.dockwidget.mFieldComboBox_dip.currentText()
        dip_dir_col = self.dockwidget.mFieldComboBox_dip_dir.currentText()
        if not dip_col or not dip_dir_col:
            self.iface.messageBar().pushMessage(
                "Please define dip & dip dir/strike fields before continuing!",
                level=Qgis.Warning, duration=15)
            return False
        return True

    def subsamplePoints(self):
        if self.dockwidget.radioButton_stochastic.isChecked():
            self.stochastic()
            return
        if not self._validate_dip_fields(require_points_layer=True):
            return
        if self.dockwidget.radioButton_gsa.isChecked():
            self.gridCellAveraging()
        elif self.dockwidget.radioButton_kent.isChecked():
            self.kent()
        elif self.dockwidget.radioButton_kent_outlier.isChecked():
            self.kentOutlier()

    def processFaults(self):
        if self.dockwidget.radioButton_fault_length.isChecked():
            self.fault_Lengths()
        elif self.dockwidget.radioButton_fault_graph.isChecked():
            self.fault_Graph()
        elif self.dockwidget.radioButton_fault_strat_offset.isChecked():
            self.fault_strat_offset()
        elif self.dockwidget.radioButton_fault_clusters.isChecked():
            self.fault_ClusterOrientations()

    def _tp_prepare_fault_attrs(self, gdf, method, fault_layer, polygon_layer,
                                strat1, strat2, strat3, strat4):
        """
        Ensure gdf has the attribute column(s) required by the chosen TP fault
        method.  Runs the relevant calculation automatically if the column is
        absent.  Returns the (possibly enriched) GeoDataFrame unchanged when
        no pre-computation is needed or when it fails.
        """
        if method == 'graph' and 'edge_type' not in gdf.columns:
            tmpdir = tempfile.mkdtemp()
            try:
                tmp_in = os.path.join(tmpdir, 'tp_graph.shp')
                gdf.to_file(tmp_in)
                FaultsGraph().CalcFaultsGraph(tmp_in)
                # save_graph_to_shapefile prepends prefix_ to the basename:
                # tp_graph_edges.shp → simplified_full_tp_graph_edges.shp
                tmp_edges = os.path.join(tmpdir, 'simplified_full_tp_graph_edges.shp')
                if os.path.exists(tmp_edges):
                    return gpd.read_file(tmp_edges).set_crs(gdf.crs, allow_override=True)
            except Exception as exc:
                print(f"T&P: graph pre-compute failed: {exc}")
            finally:
                shutil.rmtree(tmpdir, ignore_errors=True)

        elif method == 'clusters' and not any('cluster' in c.lower() for c in gdf.columns):
            tmpdir = tempfile.mkdtemp()
            try:
                tmp_in = os.path.join(tmpdir, 'tp_cl.shp')
                tmp_az = os.path.join(tmpdir, 'tp_cl_az.shp')
                gdf.to_file(tmp_in)
                fo = FaultsOrientations()
                fo.add_endpoint_azimuth(tmp_in, tmp_az, azimuth_field='endpt_az')
                best_gdf, best_score = None, -1
                for n in range(2, 8):
                    try:
                        res = fo.analyze_shapefile_trends(
                            shapefile_path=tmp_az,
                            trend_field='endpt_az',
                            n_clusters=n,
                        )
                        sc = res.get('silhouette_score', -1)
                        if sc > best_score:
                            best_score, best_gdf = sc, res.get('gdf_with_clusters')
                    except Exception:
                        pass
                if best_gdf is not None:
                    return best_gdf.set_crs(gdf.crs, allow_override=True)
            except Exception as exc:
                print(f"T&P: cluster pre-compute failed: {exc}")
            finally:
                shutil.rmtree(tmpdir, ignore_errors=True)

        elif method == 'strat_offset':
            has_strat = any('strat' in c.lower() or 'offset' in c.lower()
                            for c in gdf.columns)
            strat_columns = [c for c in [strat1, strat2, strat3, strat4] if c]
            if (not has_strat and strat_columns
                    and polygon_layer and os.path.exists(self._resolved_source(polygon_layer))):
                tmpdir = tempfile.mkdtemp()
                try:
                    tmp_in  = os.path.join(tmpdir, 'tp_so.shp')
                    tmp_out = os.path.join(tmpdir, 'tp_so_out.shp')
                    gdf.to_file(tmp_in)
                    FaultStratOffset().CalcFaultStratOffset(
                        tmp_in, self._resolved_source(polygon_layer), tmp_out,
                        strat_columns, offset_distance=50)
                    if os.path.exists(tmp_out):
                        return gpd.read_file(tmp_out).set_crs(gdf.crs, allow_override=True)
                except Exception as exc:
                    print(f"T&P: strat-offset pre-compute failed: {exc}")
                finally:
                    shutil.rmtree(tmpdir, ignore_errors=True)

        return gdf

    def topferPillewizer(self):
        """
        Apply Töpfer & Pillewizer (1966) scaling to all available layers.

        TN = ON * (OS/TS)^(x/2)   x=1 points, x=2 lines, x=3 polygons

        The OS/TS ratio and optional increment are read from the T&P group.
        If increment < ratio the method iterates through intermediate ratios
        with step sizes growing by 1.5× each time.
        """
        if not self._validate_dip_fields(require_points_layer=False):
            return

        try:
            ratio     = float(self.dockwidget.lineEdit_tp_ratio.text())
            increment = float(self.dockwidget.lineEdit_tp_increment.text())
        except ValueError:
            self.iface.messageBar().pushMessage(
                "T&P Scaling: enter valid numbers for ratio and increment.",
                level=Qgis.Warning, duration=10)
            return

        if ratio <= 0 or ratio > 1:
            self.iface.messageBar().pushMessage(
                "T&P Scaling: OS/TS ratio must be in the range (0, 1].",
                level=Qgis.Warning, duration=10)
            return

        iterations = scale_iterations(ratio, increment)

        # --- Determine selected orientation method -----------------------
        if self.dockwidget.radioButton_stochastic.isChecked():
            point_method = 'stochastic'
        elif self.dockwidget.radioButton_gsa.isChecked():
            point_method = 'gridcell_average'
        elif self.dockwidget.radioButton_kent.isChecked():
            point_method = 'spherical_kent'
        else:
            point_method = 'outlier_removal'

        # --- Determine selected fault method -----------------------------
        if self.dockwidget.radioButton_fault_length.isChecked():
            fault_method = 'length'
        elif self.dockwidget.radioButton_fault_graph.isChecked():
            fault_method = 'graph'
        elif self.dockwidget.radioButton_fault_strat_offset.isChecked():
            fault_method = 'strat_offset'
        else:
            fault_method = 'clusters'

        # --- Polygon parameters (shared with existing tools) -------------
        strat1    = self.dockwidget.mFieldComboBox_priority_1.currentText()
        strat2    = self.dockwidget.mFieldComboBox_priority_2.currentText()
        strat3    = self.dockwidget.mFieldComboBox_priority_3.currentText()
        strat4    = self.dockwidget.mFieldComboBox_priority_4.currentText()
        lithoname = self.dockwidget.mFieldComboBox_priority_5.currentText()
        dyke_field = self.dockwidget.mFieldComboBox_dyke.currentText()
        dyke_codes = (
            self.dockwidget.plainTextEdit_dyke_Codes.toPlainText()
            .replace(" ", "").split(",")
        )
        try:
            dist_thresh = float(self.dockwidget.lineEdit_node_tolerance.text())
        except ValueError:
            dist_thresh = 1.0

        fault_layer   = self.dockwidget.mMapLayerComboBox_fault_polylines.currentLayer()
        polygon_layer = self.dockwidget.mMapLayerComboBox_maps_polygons.currentLayer()

        # Field combo boxes reflect the original layer's field names, which may
        # have been shortened by _resolved_source's shapefile export (DBF field
        # names are capped at 10 characters).
        if polygon_layer:
            strat1     = self._resolved_field(polygon_layer, strat1)
            strat2     = self._resolved_field(polygon_layer, strat2)
            strat3     = self._resolved_field(polygon_layer, strat3)
            strat4     = self._resolved_field(polygon_layer, strat4)
            lithoname  = self._resolved_field(polygon_layer, lithoname)
            dyke_field = self._resolved_field(polygon_layer, dyke_field)

        # Pre-load GeoDataFrames once; each iteration feeds into the next
        current_fault_gdf   = None
        current_polygon_gdf = None

        if fault_layer and os.path.exists(self._resolved_source(fault_layer)):
            current_fault_gdf = gpd.read_file(self._resolved_source(fault_layer))
            current_fault_gdf = self._tp_prepare_fault_attrs(
                current_fault_gdf, fault_method, fault_layer, polygon_layer,
                strat1, strat2, strat3, strat4)

        if polygon_layer and os.path.exists(self._resolved_source(polygon_layer)):
            current_polygon_gdf = gpd.read_file(self._resolved_source(polygon_layer))
            crs = polygon_layer.crs()
            if crs.isGeographic():
                dist_thresh = dist_thresh / 110000

        # Points: set up engine and columns once; chain the GDF across iterations
        current_pt_gdf  = None
        pt_engine       = None
        pt_dip_col      = None
        pt_dip_dir_col  = None
        pt_result = self.setUpPointSampler()
        if pt_result:
            current_pt_gdf, pt_engine, pt_dip_col, pt_dip_dir_col = pt_result

        try:
            grid_sz = float(self.dockwidget.lineEdit_grid_size.text())
        except ValueError:
            grid_sz = 5000.0

        # Points are always subsampled from the original full dataset using the
        # absolute cumulative ratio, so results at each scale are independent.
        # Lines and polygons chain incrementally (each step feeds the next).
        original_pt_gdf = current_pt_gdf  # never mutated

        # --- Ratio=1 baseline: all layers at full extent ---
        if original_pt_gdf is not None:
            self.finalisePointSampler(original_pt_gdf, self.points_layer, "tp_1")

        if current_fault_gdf is not None:
            ln_full = scale_lines_tp(current_fault_gdf, 1.0, fault_method)
            layer_path = os.path.dirname(self._resolved_source(fault_layer))
            out_name = f"{fault_layer.name()}_tp_1.shp"
            out_path = os.path.join(layer_path, out_name)
            if os.path.exists(out_path):
                out_path = out_path.replace(
                    '.shp', f"_{random.randint(10000,99999)}.shp")
            ln_full.to_file(out_path, driver='ESRI Shapefile')
            ln_layer = QgsVectorLayer(out_path, f"{fault_layer.name()}_tp_1", "ogr")
            self._add_layer_with_style(ln_layer, fault_layer)

        if current_polygon_gdf is not None:
            layer_path = os.path.dirname(self._resolved_source(polygon_layer))
            out_name = f"{polygon_layer.name()}_tp_1.shp"
            out_path = os.path.join(layer_path, out_name)
            if os.path.exists(out_path):
                out_path = out_path.replace(
                    '.shp', f"_{random.randint(10000,99999)}.shp")
            current_polygon_gdf.to_file(out_path, driver='ESRI Shapefile')
            poly_layer = QgsVectorLayer(out_path, f"{polygon_layer.name()}_tp_1", "ogr")
            self._add_layer_with_style(poly_layer, polygon_layer)

        prev_ratio = 1.0
        for iter_ratio in iterations:
            step_ratio = iter_ratio / prev_ratio   # still needed for lines/polygons
            prev_ratio = iter_ratio
            ratio_tag  = f"{iter_ratio:.6f}".rstrip('0').rstrip('.')

            if abs(iter_ratio - 1.0) < 1e-9:
                continue  # ratio=1 baseline already output above

            # --- Points (always from original, absolute ratio) ---
            if original_pt_gdf is not None:
                pt_scaled = scale_points_tp(
                    original_pt_gdf, iter_ratio, pt_engine, point_method,
                    pt_dip_col, pt_dip_dir_col, grid_sz)

                # Drop internal helper columns for the saved file
                drop = [c for c in ['EASTING', 'NORTHING', '_eff_dipdir']
                        if c in pt_scaled.columns]
                if drop:
                    pt_scaled = pt_scaled.drop(columns=drop)

                self.finalisePointSampler(
                    pt_scaled, self.points_layer, f"tp_{ratio_tag}")

            # --- Lines ---
            if current_fault_gdf is not None:
                if fault_method == 'graph':
                    # Removing faults changes node topology, so edge_type must be
                    # recomputed from the surviving fault set each iteration.
                    if 'edge_type' in current_fault_gdf.columns:
                        current_fault_gdf = current_fault_gdf.drop(columns=['edge_type'])
                    current_fault_gdf = self._tp_prepare_fault_attrs(
                        current_fault_gdf, fault_method, fault_layer, polygon_layer,
                        strat1, strat2, strat3, strat4)
                ln_scaled = scale_lines_tp(current_fault_gdf, step_ratio, fault_method)
                current_fault_gdf = ln_scaled
                layer_path = os.path.dirname(self._resolved_source(fault_layer))
                out_name   = f"{fault_layer.name()}_tp_{ratio_tag}.shp"
                out_path   = os.path.join(layer_path, out_name)
                if os.path.exists(out_path):
                    out_path = out_path.replace(
                        '.shp', f"_{random.randint(10000,99999)}.shp")
                ln_scaled.to_file(out_path, driver='ESRI Shapefile')
                ln_layer = QgsVectorLayer(
                    out_path,
                    f"{fault_layer.name()}_tp_{ratio_tag}", "ogr")
                self._add_layer_with_style(ln_layer, fault_layer)

            # --- Polygons ---
            if current_polygon_gdf is not None:
                try:
                    poly_scaled = scale_polygons_tp(
                        current_polygon_gdf, step_ratio,
                        StructuralPolygonSubSampler,
                        distance_threshold=dist_thresh,
                        lithoname=lithoname or None,
                        strat1=strat1 or None,
                        strat2=strat2 or None,
                        strat3=strat3 or None,
                        strat4=strat4 or None,
                        dyke_field=dyke_field or None,
                        dyke_codes=dyke_codes if dyke_field else None,
                        triangulator_class=(
                            PolygonTriangulator if dyke_field else None),
                    )
                    current_polygon_gdf = poly_scaled
                except Exception as exc:
                    print(f"T&P polygon scaling error at ratio {iter_ratio}: {exc}")
                    poly_scaled = current_polygon_gdf.copy()

                layer_path = os.path.dirname(self._resolved_source(polygon_layer))
                out_name   = f"{polygon_layer.name()}_tp_{ratio_tag}.shp"
                out_path   = os.path.join(layer_path, out_name)
                if os.path.exists(out_path):
                    out_path = out_path.replace(
                        '.shp', f"_{random.randint(10000,99999)}.shp")
                poly_scaled.to_file(out_path, driver='ESRI Shapefile')
                poly_layer = QgsVectorLayer(
                    out_path,
                    f"{polygon_layer.name()}_tp_{ratio_tag}", "ogr")
                self._add_layer_with_style(poly_layer, polygon_layer)

    def firstOrder(self):
        self.polygon_layer = (
            self.dockwidget.mMapLayerComboBox_maps_polygons.currentLayer()
        )
        self.point_layer = self.dockwidget.mMapLayerComboBox_points.currentLayer()
        if self.polygon_layer is not None and self.point_layer is not None:
            if os.path.exists(self._resolved_source(self.polygon_layer)) and os.path.exists(self._resolved_source(self.point_layer)):
                distance_threshold = float(
                    self.dockwidget.lineEdit_1o_distance.text()
                )
                angle_threshold = float(self.dockwidget.lineEdit_1o_angle.text())

                # Convert QGIS layer to GeoPandas
                geology         = gpd.read_file(self._resolved_source(self.polygon_layer))
                contact_gdf     = geology[['geometry']].copy()
                contact_gdf['geometry'] = geology.boundary
                orientation_gdf = gpd.read_file(self._resolved_source(self.point_layer))
                dip_col = self.dockwidget.mFieldComboBox_dip.currentText()
                dip_dir_col = self.dockwidget.mFieldComboBox_dip_dir.currentText()
                if not dip_col or not dip_dir_col:
                    print("First Order: please select DIP and DIP_DIR fields before running.")
                    return
                dip_col = self._resolved_field(self.point_layer, dip_col)
                dip_dir_col = self._resolved_field(self.point_layer, dip_dir_col)

                # Create an instance of SubsamplingEngine
                subsampler = SubsamplingEngine(dip=dip_col, dipdir=dip_dir_col)
                kwargs={}
                kwargs['contact_gdf'] = contact_gdf
                kwargs['dist_buffer']= distance_threshold
                kwargs['angle_tol']= angle_threshold


                output_gdf=subsampler.subsample("firstorder", orientation_gdf, path_out="", **kwargs)


                # Write GeoPandas back to file and reload
                self.finalisePointSampler(
                    output_gdf, self.point_layer, f"first_order_{int(distance_threshold)}_{int(angle_threshold)}")
    def minPolyArea(self):
        self.polygon_layer = (
            self.dockwidget.mMapLayerComboBox_maps_polygons.currentLayer()
        )
        if self.polygon_layer is not None:
            if os.path.exists(self._resolved_source(self.polygon_layer)):
                distance_threshold = float(
                    self.dockwidget.lineEdit_node_tolerance.text()
                )
                lithoname = self.dockwidget.mFieldComboBox_priority_5.currentText()
                strat1 = self.dockwidget.mFieldComboBox_priority_1.currentText()
                strat2 = self.dockwidget.mFieldComboBox_priority_2.currentText()
                strat3 = self.dockwidget.mFieldComboBox_priority_3.currentText()
                strat4 = self.dockwidget.mFieldComboBox_priority_4.currentText()
                dyke_codes = (
                    self.dockwidget.plainTextEdit_dyke_Codes.toPlainText()
                    .replace(" ", "")
                    .split(",")
                )
                dyke_field = self.dockwidget.mFieldComboBox_dyke.currentText()
                dyke_index = self.dockwidget.mFieldComboBox_dyke.currentIndex()
                lithoname  = self._resolved_field(self.polygon_layer, lithoname)
                strat1     = self._resolved_field(self.polygon_layer, strat1)
                strat2     = self._resolved_field(self.polygon_layer, strat2)
                strat3     = self._resolved_field(self.polygon_layer, strat3)
                strat4     = self._resolved_field(self.polygon_layer, strat4)
                dyke_field = self._resolved_field(self.polygon_layer, dyke_field)
                incScale = self.dockwidget.mQgsDoubleSpinBox_upinc.value()

                maxScale = float(self.dockwidget.lineEdit_polygon_area.text())
                if self.dockwidget.checkBox_series.isChecked():
                    minScale = 0.0

                else:
                    incScale = maxScale 
                    minScale = maxScale

                for upScale in np.arange(minScale, maxScale + 0.001, incScale):
                    min_area_threshold = ((upScale / 2) ** 2) * np.pi
                    parameter = str(upScale)
                    crs = self.polygon_layer.crs()
                    if crs.isGeographic():
                        min_area_threshold = min_area_threshold / (110000 * 110000)
                        distance_threshold = distance_threshold / 110000

                    layer_path = os.path.dirname(self._resolved_source(self.polygon_layer))
                    new_path = (
                        layer_path
                        + "/"
                        + self.polygon_layer.name()
                        + "_min_area_"
                        + parameter
                        + ".shp"
                    )
                    if os.path.exists(new_path):
                        random_5_digit_integer = random.randint(10000, 99999)
                        new_path = (
                            layer_path
                            + "/"
                            + self.polygon_layer.name()
                            + "_min_area_"
                            + parameter
                            + "_"
                            + str(random_5_digit_integer)
                            + ".shp"
                        )

                    # Convert QGIS layer to GeoPandas if first in series
                    if upScale == minScale:
                        gdf = gpd.read_file(self._resolved_source(self.polygon_layer))

                    # Handle dykes special case
                    if dyke_index != 0:
                        triangulator = PolygonTriangulator(
                            gdf=gdf,
                            id_column=dyke_field,
                            min_area_threshold=min_area_threshold * 1000000.0,
                            distance_threshold=distance_threshold,
                            strat1=strat1,
                            strat2=strat2,
                            strat3=strat3,
                            strat4=strat4,
                            lithoname=lithoname,
                        )
                        gdf = triangulator.triangulate_polygons(target_ids=dyke_codes)

                    # Create an instance of StructuralPolygonSubSampler
                    polygonSubsampler = StructuralPolygonSubSampler(gdf)
                    # Modified parameter order: stratigraphic fields come before lithological field
                    output_gdf = polygonSubsampler.clean_small_polygons_and_holes_new(
                        gdf,
                        min_area_threshold=min_area_threshold * 1000000.0,
                        distance_threshold=distance_threshold,
                        strat1=strat1,
                        strat2=strat2,
                        strat3=strat3,
                        strat4=strat4,
                        lithoname=lithoname,
                    )
                    output_gdf.to_file(new_path, driver="ESRI Shapefile")

                    # reload as layer
                    upscaled_layer = QgsVectorLayer(
                        new_path,
                        self.polygon_layer.name() + "_min_area_" + parameter,
                        "ogr",
                    )

                    self._add_layer_with_style(upscaled_layer, self.polygon_layer,
                                               self.polygon_layer.name() + "_min_area")
            else:
                self.iface.messageBar().pushMessage(
                    "Sorry, only layers saved to disk as a shapefile can be processed!",
                    level=Qgis.Warning,
                    duration=15,
                )

    def mergeSegments(self):
        self.polyline_layer = (
            self.dockwidget.mMapLayerComboBox_fault_polylines.currentLayer()
        )
        if os.path.exists(self._resolved_source(self.polyline_layer)):
            distance_tolerance = float(self.dockwidget.lineEdit_merge_tolerance.text())
            angle_tolerance = float(self.dockwidget.lineEdit_merge_search_angle.text())
            min_join_angle = float(self.dockwidget.lineEdit_merge_join_angle.text())

            merger = FaultLineMerger(
                distance_tolerance=distance_tolerance,  # 10 map units
                angle_tolerance=angle_tolerance,  # 30 degrees for line parallelism
                min_join_angle=min_join_angle,  # Minimum 150° angle at join point (rejects angles < 150°)
            )

            layer_path = os.path.dirname(self._resolved_source(self.polyline_layer))
            new_path = layer_path + "/" + self.polyline_layer.name() + "_merged.shp"
            if os.path.exists(new_path):
                random_5_digit_integer = random.randint(10000, 99999)
                new_path = (
                    layer_path
                    + "/"
                    + self.polyline_layer.name()
                    + "_merged_"
                    + str(random_5_digit_integer)
                    + ".shp"
                )

            gdf = gpd.read_file(self._resolved_source(self.polyline_layer))
            output_gdf = merger.process_shapefile(gdf)
            output_gdf.to_file(new_path, driver="ESRI Shapefile")

            # reload as layer
            upscaled_layer = QgsVectorLayer(
                new_path,
                self.polyline_layer.name() + "_merged",
                "ogr",
            )

            self._add_layer_with_style(upscaled_layer, self.polyline_layer,
                                       self.polyline_layer.name() + "_merge")

        else:
            self.iface.messageBar().pushMessage(
                "Sorry, only layers saved to disk as a shapefile can be processed!",
                level=Qgis.Warning,
                duration=15,
            )

    def fault_Graph(self):
        self.polyline_layer = (
            self.dockwidget.mMapLayerComboBox_fault_polylines.currentLayer()
        )
        if os.path.exists(self._resolved_source(self.polyline_layer)):
            graph_obj = FaultsGraph()
            graph_obj.CalcFaultsGraph(self._resolved_source(self.polyline_layer))
            directory, filename = os.path.split(self._resolved_source(self.polyline_layer))
            new_path = (
                directory
                + "/simplified_full_"
                + str(filename).replace(".shp", "_edges.shp")
            )
            print("dir fn np", directory, filename, new_path)
            graph_layer = QgsVectorLayer(
                new_path, self.polyline_layer.name() + "_graph", "ogr"
            )

            self._add_layer_with_style(graph_layer, self.polyline_layer,
                                       self.polyline_layer.name() + "_graph")

    def fault_strat_offset(self):
        self.polyline_layer = (
            self.dockwidget.mMapLayerComboBox_fault_polylines.currentLayer()
        )
        if os.path.exists(self._resolved_source(self.polyline_layer)) and os.path.exists(
            self._resolved_source(self.map_layer)
        ):

            strat_columns = []
            if self.dockwidget.mFieldComboBox_priority_1.currentText():
                strat_columns.append(
                    self.dockwidget.mFieldComboBox_priority_1.currentText()
                )
            if self.dockwidget.mFieldComboBox_priority_2.currentText():
                strat_columns.append(
                    self.dockwidget.mFieldComboBox_priority_2.currentText()
                )
            if self.dockwidget.mFieldComboBox_priority_3.currentText():
                strat_columns.append(
                    self.dockwidget.mFieldComboBox_priority_3.currentText()
                )
            if self.dockwidget.mFieldComboBox_priority_4.currentText():
                strat_columns.append(
                    self.dockwidget.mFieldComboBox_priority_4.currentText()
                )

            if not strat_columns:
                print("Fault Strat Offset: please select at least one priority field before running.")
                return

            strat_columns = [self._resolved_field(self.map_layer, c) for c in strat_columns]

            strat_offset_obj = FaultStratOffset()
            strat_offset_obj.CalcFaultStratOffset(
                self._resolved_source(self.polyline_layer),
                self._resolved_source(self.map_layer),
                self._resolved_source(self.polyline_layer).replace(".shp", "_stratOffset.shp"),
                strat_columns,
                offset_distance=50,
            )

            directory, filename = os.path.split(self._resolved_source(self.polyline_layer))
            new_path = (
                directory + "/" + str(filename).replace(".shp", "_stratOffset.shp")
            )
            print("dir fn np", directory, filename, new_path)
            graph_layer = QgsVectorLayer(
                new_path, self.polyline_layer.name() + "_stratOffset", "ogr"
            )

            self._add_layer_with_style(graph_layer, self.polyline_layer,
                                       self.polyline_layer.name() + "_stratOffset")

    def fault_Lengths(self):
        self.polyline_layer = (
            self.dockwidget.mMapLayerComboBox_fault_polylines.currentLayer()
        )
        if os.path.exists(self._resolved_source(self.polyline_layer)):

            crs = self.polyline_layer.crs()
            if crs.isGeographic():
                crs_units = "degrees"
            else:
                crs_units = "meters"

            outPath = self._resolved_source(self.polyline_layer).replace(".shp", "_length.shp")
            filter_obj = FaultLengths()
            gdf = filter_obj.add_polyline_length(
                self._resolved_source(self.polyline_layer),
                output_path=outPath,
                length_field="line_len",
                unit=crs_units,
                force_projection=None,
            )

            # reload as layer
            upscaled_layer = QgsVectorLayer(
                outPath,
                self.polyline_layer.name() + "_Length",
                "ogr",
            )

            self._add_layer_with_style(upscaled_layer, self.polyline_layer,
                                       self.polyline_layer.name() + "_Length")

        else:
            self.iface.messageBar().pushMessage(
                "Sorry, only layers saved to disk as a shapefile can be processed!",
                level=Qgis.Warning,
                duration=15,
            )

    def fault_ClusterOrientations(self):
        print("fault_ClusterOrientations")
        print(
            "polyline_layer:",
            self.dockwidget.mMapLayerComboBox_fault_polylines.currentText(),
        )
        self.polyline_layer = (
            self.dockwidget.mMapLayerComboBox_fault_polylines.currentLayer()
        )
        if os.path.exists(self._resolved_source(self.polyline_layer)):
            outPath = self._resolved_source(self.polyline_layer).replace(".shp", "_endpt_az.shp")
            filter_obj = FaultsOrientations()
            gdf = filter_obj.add_endpoint_azimuth(
                self._resolved_source(self.polyline_layer), outPath, azimuth_field="endpt_az"
            )

            best_n = filter_obj.example_manual_clusters(
                outPath,
                self.dockwidget.mMapLayerComboBox_fault_polylines.currentText(),
                azimuth_field="endpt_az",
            )
            layer_path = os.path.dirname(self._resolved_source(self.polyline_layer))
            new_path = (
                layer_path
                + "/"
                + self.polyline_layer.name()
                + f"_fault_clusters_{best_n}.shp"
            )

            upscaled_layer = QgsVectorLayer(
                new_path,
                self.polyline_layer.name() + f"_fault_clusters_{best_n}",
                "ogr",
            )

            self._add_layer_with_style(upscaled_layer, self.polyline_layer,
                                       self.polyline_layer.name() + f"_fault_clusters_{best_n}")

    def updatePointsFields(self):
        """Update the fields in the points layer combo box when a new layer is selected."""
        layerName = self.dockwidget.mMapLayerComboBox_points.currentText()
        if layerName:
            self.points_layer = QgsProject.instance().mapLayersByName(layerName)[0]
            self.dockwidget.mFieldComboBox_dip.setLayer(self.points_layer)
            self.dockwidget.mFieldComboBox_dip_dir.setLayer(self.points_layer)

    def simplifyMap(self):
        polygon_layer = self.dockwidget.mMapLayerComboBox_maps_polygons.currentLayer()
        fault_layer = self.dockwidget.mMapLayerComboBox_fault_polylines.currentLayer()
        if polygon_layer is None or fault_layer is None:
            return
        area_tolerance = int(self.dockwidget.lineEdit_area_tolerance.text())
        input_file = self._resolved_source(polygon_layer)
        fault_file = self._resolved_source(fault_layer)
        output_file = input_file.replace(".shp", f"_MVW_{area_tolerance}.shp")
        poly_name=polygon_layer.name()+f"_MVW_{area_tolerance}"
        fault_path=self._resolved_source(fault_layer).replace(".shp", f"_sn_si_{area_tolerance}.shp")
        fault_name=fault_layer.name()+f"_sn_si_{area_tolerance}"
        result = vector_simplify_file_two_stage(
                        input_file=input_file,
                        output_file=output_file,
                        method="modified_visvalingam_whyatt",
                        threshold=area_tolerance,
                        fault_file=fault_file
                    )
        upscaled_layer = QgsVectorLayer(
                output_file,
                poly_name,
                "ogr",
            )

        self._add_layer_with_style(upscaled_layer, polygon_layer, poly_name)

        upscaled_fault = QgsVectorLayer(fault_path, fault_name, "ogr")
        self._add_layer_with_style(upscaled_fault, fault_layer, fault_name)

    def updateMapsFields(self):
        """Update the fields in the maps layer combo box when a new layer is selected."""
        layerName = self.dockwidget.mMapLayerComboBox_maps_polygons.currentText()
        if layerName:
            self.map_layer = QgsProject.instance().mapLayersByName(layerName)[0]
            self.dockwidget.mFieldComboBox_priority_1.setLayer(self.map_layer)
            self.dockwidget.mFieldComboBox_priority_2.setLayer(self.map_layer)
            self.dockwidget.mFieldComboBox_priority_3.setLayer(self.map_layer)
            self.dockwidget.mFieldComboBox_priority_4.setLayer(self.map_layer)
            self.dockwidget.mFieldComboBox_priority_5.setLayer(self.map_layer)
            self.dockwidget.mFieldComboBox_dyke.setLayer(self.map_layer)
            self.dockwidget.mFieldComboBox_priority_1.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_priority_2.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_priority_3.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_priority_4.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_priority_5.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_dyke.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_priority_1.setCurrentIndex(0)
            self.dockwidget.mFieldComboBox_priority_2.setCurrentIndex(0)
            self.dockwidget.mFieldComboBox_priority_3.setCurrentIndex(0)
            self.dockwidget.mFieldComboBox_priority_4.setCurrentIndex(0)
            self.dockwidget.mFieldComboBox_priority_5.setCurrentIndex(0)
            self.dockwidget.mFieldComboBox_dyke.setCurrentIndex(0)

    def defineTips(self):
        """
        Defines tooltips for various UI elements in the QGIS plugin.
        This method assigns descriptive tooltips to different widgets and controls
        in the plugin's user interface. These tooltips provide guidance and
        explanations for the functionality of each element, helping users
        understand their purpose and usage.
        """

        self.dockwidget.mMapLayerComboBox_points.setToolTip(
            "Points layer selected for processing"
        )
        self.dockwidget.mFieldComboBox_dip.setToolTip(
            "Field containing Dip information"
        )
        self.dockwidget.mFieldComboBox_dip_dir.setToolTip(
            "Field containing Dip Direction or Strike information"
        )
        self.dockwidget.checkBox_dip_dir.setToolTip(
            "Dip Direction used instead of Strike \nProbably not needed?"
        )
        self.dockwidget.radioButton_stochastic.setToolTip(
            "Keep a randomly selected percentage of points"
        )
        self.dockwidget.mQgsDoubleSpinBox_percent.setToolTip(
            "Percentage of points to keep, between 0 and 100"
        )
        self.dockwidget.radioButton_gsa.setToolTip(
            "Grid Cell Averaging of points"
        )
        self.dockwidget.lineEdit_grid_size.setToolTip(
            "Grid size in layer units for cell averaging"
        )
        self.dockwidget.radioButton_kent.setToolTip(
            "Grid Cell Averaging of points using Kent statistics"
        )
        self.dockwidget.lineEdit_grid_size_kent.setToolTip(
            "Grid size in layer units for cell averaging using Kent statistics"
        )
        self.dockwidget.radioButton_kent_outlier.setToolTip(
            "Grid Cell Averaging using Kent statistics, removing biggest outlier"
        )
        self.dockwidget.lineEdit_grid_size_kent_2.setToolTip(
            "Grid size in layer units for Kent outlier averaging"
        )
        self.dockwidget.pushButton_subsample_points.setToolTip(
            "Run the selected subsampling method on the point layer"
        )

        self.dockwidget.pushButton_1o_sampling.setToolTip(
            "First order retention of points based on distance and angle to contacts"
        )
        self.dockwidget.lineEdit_1o_distance.setToolTip(
            "Distance buffer for first order retention of points"
        )
        self.dockwidget.lineEdit_1o_angle.setToolTip(
            "Angle tolerance for first order retention of points"
        )
        self.dockwidget.mMapLayerComboBox_maps_polygons.setToolTip(
            "Layer map polygons selected for processing"
        )
        self.dockwidget.lineEdit_node_tolerance.setToolTip(
            "Tolerance for node snapping in layer units"
        )
        self.dockwidget.mFieldComboBox_priority_1.setToolTip(
            "Field selected for 1st priority of polygon merging and fault strat offset"
        )
        self.dockwidget.mFieldComboBox_priority_2.setToolTip(
            "Field selected for 2nd priority of polygon merging and fault strat offset"
        )
        self.dockwidget.mFieldComboBox_priority_3.setToolTip(
            "Field selected for 3rd priority of polygon merging and fault strat offset"
        )
        self.dockwidget.mFieldComboBox_priority_4.setToolTip(
            "Field selected for 4th priority of polygon merging and fault strat offset"
        )
        self.dockwidget.mFieldComboBox_priority_5.setToolTip(
            "Field selected for 5th priority of polygon merging and fault strat offset"
        )
        self.dockwidget.mFieldComboBox_dyke.setToolTip("Field selected for dyke codes")
        self.dockwidget.plainTextEdit_dyke_Codes.setToolTip(
            "Comma separated list of dyke codes for special handling"
        )
        self.dockwidget.lineEdit_polygon_area.setToolTip(
            "Minimum threshold diameter of polygons to be retained (in km)"
        )
        self.dockwidget.checkBox_series.setToolTip(
            "Calculate full series of rescaling maps between increment and threshold"
        )
        self.dockwidget.mQgsDoubleSpinBox_upinc.setToolTip(
            "Increment for arithmetic series of maps"
        )
        self.dockwidget.pushButton_minPolyArea.setToolTip("Rescale map polygons")
        self.dockwidget.mMapLayerComboBox_fault_polylines.setToolTip(
            "Fault polyline layer selected for processing"
        )
        self.dockwidget.pushButton_merge_segments.setToolTip(
            "Merge connected polylines\nWARNING: VERY SLOW FOR LARGE DATASETS (>1s per object)"
        )
        self.dockwidget.lineEdit_merge_tolerance.setToolTip(
            "Minimum distance between end points"
        )
        self.dockwidget.lineEdit_merge_search_angle.setToolTip(
            "Minimum angle mismatch between end segments"
        )
        self.dockwidget.lineEdit_merge_join_angle.setToolTip(
            "Maximum angle for newly-formed angle defined by end segments"
        )
        self.dockwidget.radioButton_fault_length.setToolTip(
            "Add fault length attribute\n"
            "Does not filter data, just adds length attribute"
        )
        self.dockwidget.radioButton_fault_graph.setToolTip(
            "Add graph data to faults\n"
            "Does not filter data, just adds graph attributes"
        )
        self.dockwidget.radioButton_fault_clusters.setToolTip(
            "Add suite of orientation cluster data to faults (2-7 clusters)\n"
            "Does not filter data, just adds cluster attributes"
        )
        self.dockwidget.radioButton_fault_strat_offset.setToolTip(
            "Add stratigraphic offset data to faults\n"
            "Requires both fault and map polygon layers\n"
            "as well as Geology Map Priority Codes (see below) to be defined\n"
            "Does not filter data, just adds Strat Offset attributes"
        )
        self.dockwidget.pushButton_process_faults.setToolTip(
            "Run the selected fault attribute method"
        )
        self.dockwidget.lineEdit_tp_ratio.setToolTip(
            "OS/TS ratio — e.g. 0.1 for going from 1:50 000 to 1:500 000"
        )
        self.dockwidget.lineEdit_tp_increment.setToolTip(
            "Fractional reduction per step (0–1).\n"
            "Step n ratio = (1 - increment)^n\n"
            "e.g. increment=0.1 → steps 0.9, 0.81, 0.729 …\n"
            "Stops after the first step that drops below the OS/TS ratio."
        )
        self.dockwidget.pushButton_tp_run.setToolTip(
            "Apply Töpfer & Pillewizer scaling to all selected layers.\n"
            "Uses the orientation and fault methods selected above."
        )
        self.dockwidget.lineEdit_area_tolerance.setToolTip(
            "Tolerance for area-based filtering, in m2"
        )
        self.dockwidget.pushButton_simplifyMap.setToolTip(
            "Simplify map polygons based on area tolerance, respecting fault topology"
        )
    def run(self):
        """Run method that loads and starts the plugin"""

        if not self.pluginIsActive:
            self.pluginIsActive = True

            # print ("** STARTING GeoSubSampler")

            # dockwidget may not exist if:
            #    first run of plugin
            #    removed on close (see self.onClosePlugin method)
            if self.dockwidget == None:
                # Create the dockwidget (after translation) and keep reference
                self.dockwidget = GeoSubSamplerDockWidget()

            # connect to provide cleanup on closing of dockwidget
            self.dockwidget.closingPlugin.connect(self.onClosePlugin)

            # show the dockwidget
            # TODO: fix to allow choice of dock location
            self.iface.addDockWidget(RightDockWidgetArea, self.dockwidget)
            # Find existing dock widgets in the right area
            right_docks = [
                d
                for d in self.iface.mainWindow().findChildren(QDockWidget)
                if self.iface.mainWindow().dockWidgetArea(d) == RightDockWidgetArea
            ]
            # If there are other dock widgets, tab this one with the first one found
            if right_docks:
                for dock in right_docks:
                    if dock != self.dockwidget:
                        self.iface.mainWindow().tabifyDockWidget(dock, self.dockwidget)
                        # Optionally, bring your plugin tab to the front
                        self.dockwidget.raise_()
                        break
            # Raise the docked widget above others
            self.defineTips()
            self.dockwidget.show()
            self.dockwidget.mMapLayerComboBox_points.setFilters(
                QgsMapLayerProxyModel.PointLayer
            )

            self.dockwidget.mMapLayerComboBox_maps_polygons.setFilters(
                QgsMapLayerProxyModel.PolygonLayer
            )

            self.dockwidget.mMapLayerComboBox_fault_polylines.setFilters(
                QgsMapLayerProxyModel.LineLayer
            )
            layerName = self.dockwidget.mMapLayerComboBox_points.currentText()
            if layerName:
                self.points_layer = QgsProject.instance().mapLayersByName(layerName)[0]
                self.dockwidget.mFieldComboBox_dip.setLayer(self.points_layer)
                self.dockwidget.mFieldComboBox_dip_dir.setLayer(self.points_layer)

            layerName = self.dockwidget.mMapLayerComboBox_maps_polygons.currentText()
            if layerName:
                self.updateMapsFields()

            self.dockwidget.mMapLayerComboBox_points.layerChanged.connect(
                self.updatePointsFields
            )
            self.dockwidget.mMapLayerComboBox_maps_polygons.layerChanged.connect(
                self.updateMapsFields
            )

            self.dockwidget.pushButton_subsample_points.clicked.connect(self.subsamplePoints)
            self.dockwidget.pushButton_1o_sampling.clicked.connect(self.firstOrder)

            self.dockwidget.pushButton_minPolyArea.clicked.connect(self.minPolyArea)

            self.dockwidget.pushButton_merge_segments.clicked.connect(
                self.mergeSegments
            )
            self.dockwidget.pushButton_process_faults.clicked.connect(self.processFaults)
            self.dockwidget.pushButton_tp_run.clicked.connect(self.topferPillewizer)
            self.dockwidget.pushButton_simplifyMap.clicked.connect(
                self.simplifyMap
            )
            self.dockwidget.mFieldComboBox_dip.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_dip.setCurrentIndex(0)
            self.dockwidget.mFieldComboBox_dip_dir.setAllowEmptyFieldName(True)
            self.dockwidget.mFieldComboBox_dip_dir.setCurrentIndex(0)