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 (QGIS 3.x / 4.x enum names)
from .qt_compat import (RIGHT_DOCK_AREA as RightDockWidgetArea, MSG_INFO, MSG_WARNING,
                        MSG_CRITICAL, MSG_SUCCESS, POINT_LAYER, LINE_LAYER, POLYGON_LAYER)
QT6 = hasattr(Qt, "AlignmentFlag")


# 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, dip_field_problem
from .calcs.TopferPillewizer import dipdir_to_strike, 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 csv
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", "en", type=str) or "en")[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'):
            if not os.path.exists(path_part) or self._dir_writable(os.path.dirname(path_part)):
                return path_part
            # Read-only folder: tools write results beside the source, so work on a copy
            return self._copy_to_writable(path_part)

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

        out_dir = os.path.dirname(path_part)
        if not self._dir_writable(out_dir):
            out_dir = self._fallback_dir()
        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))
            fields_selected = True

            if (dip_col == "" or dip_dir_col == "") and checkFields:
                self.iface.messageBar().pushMessage(
                    "Please define dip & dip dir/strike fields before continuing!",
                    level=MSG_WARNING,
                    duration=15,
                )
                return False
            elif dip_col == "" or dip_dir_col == "":
                fields_selected = False
                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=MSG_CRITICAL, duration=20)
                return False

            # Dip must be 0-90: values above that mean the dip and dip dir/strike fields
            # have been swapped (only checked for fields the user actually selected)
            if fields_selected:
                problem = dip_field_problem(gdf[dip_col], dip_col)
                if problem:
                    self._warn(problem + " Nothing was processed.", level=MSG_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=MSG_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=MSG_WARNING,
                duration=15,
            )
            return False

    @staticmethod
    def _parse_codes(text):
        """Split a comma/newline separated code list, trimming whitespace and quotes."""
        codes = []
        for part in text.replace("\n", ",").replace(";", ",").split(","):
            part = part.strip().strip("'\"").strip()
            if part:
                codes.append(part)
        return codes

    def _report_ignore_codes(self, triangulator, codes, gdf):
        """Tell the user, at the top of the pane, what the ignore codes matched."""
        s = getattr(triangulator, 'last_summary', None)
        if s is None:
            return
        if s['codes_not_found']:
            available = sorted(set(str(v) for v in gdf[triangulator.id_column].unique()))
            shown = ", ".join(available[:15]) + (" ..." if len(available) > 15 else "")
            self._warn(
                f"Ignore codes not found in field '{triangulator.id_column}': "
                f"{', '.join(s['codes_not_found'])}. Values present: {shown}")
        if s.get('failed'):
            self._warn(
                f"Ignore codes: {s['failed']} polygon(s) could not be triangulated and were "
                f"left unchanged (see console).")
        if s['matched'] and not s['triangulated'] and not s.get('failed'):
            self._warn(
                f"Ignore codes matched {s['matched']} polygon(s), but none is smaller than "
                f"the area threshold ({s['threshold']:.4g} map units²), so none was "
                f"triangulated. Only polygons below the threshold are triangulated.",
                duration=15)
        elif s['triangulated']:
            self.iface.messageBar().pushMessage(
                f"Ignore codes: {s['triangulated']} of {s['matched']} matching polygon(s) "
                f"triangulated ({s['above_threshold']} above the area threshold left as is).",
                level=MSG_INFO, duration=8)

    def _warn(self, message, level=MSG_WARNING, duration=15):
        """Print to the console and show the message at the top of the QGIS pane."""
        print(message)
        self.iface.messageBar().pushMessage(message, level=level, duration=duration)

    def _guard(self, fn, label):
        """
        Wrap a button handler so an unexpected exception is reported at the top
        of the pane (and traceback printed to the console) instead of failing
        silently. Qt's clicked(bool) argument is deliberately dropped.
        """
        def wrapper(*_args):
            try:
                fn()
            except Exception as exc:
                import traceback
                traceback.print_exc()
                detail = str(exc)
                if "NoneType" in detail:
                    detail = "a required layer is not selected (" + detail + ")"
                self._warn(f"{label} failed: {detail}", level=MSG_CRITICAL, duration=20)
        return wrapper

    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():
            self._warn(
                f"Failed to load result layer '{fail_msg or new_layer.name()}' - "
                f"the tool produced no output (see console).",
                level=MSG_CRITICAL)
            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)

    @staticmethod
    def _safe_filename(text, max_len=80):
        """
        Make a layer name safe to use inside a file name. Characters that are
        illegal in Windows file names (\\ / : * ? " < > |) and control characters
        make the write fail even in a writable folder, e.g. a layer called
        "Structure: 2024" or "points/clean".
        """
        cleaned = re.sub(r'[\\/:*?"<>|\x00-\x1f]+', '_', str(text))
        cleaned = cleaned.strip(' .')
        return cleaned[:max_len] or 'layer'

    def _output_directory(self, qgis_layer):
        """
        Folder for results: next to the layer's file when that exists and can be
        written to, otherwise the project folder, otherwise the temp folder.
        Memory/temporary layers and web sources have no usable folder.
        """
        resolved = self._resolved_source(qgis_layer)
        folder = os.path.dirname(resolved.split('|')[0])
        return folder if self._dir_writable(folder) else self._fallback_dir()

    def _dir_writable(self, folder):
        """True if a file can really be created in `folder` (os.access is unreliable on Windows)."""
        if not folder or not os.path.isdir(folder):
            return False
        cache = self.__dict__.setdefault('_writable_cache', {})
        if folder not in cache:
            try:
                with tempfile.NamedTemporaryFile(dir=folder):
                    pass
                cache[folder] = True
            except OSError:
                cache[folder] = False
        return cache[folder]

    def _fallback_dir(self):
        """Project folder if writable, otherwise the temp folder."""
        for folder in (QgsProject.instance().homePath(), tempfile.gettempdir()):
            if self._dir_writable(folder):
                return folder
        return tempfile.gettempdir()

    def _copy_to_writable(self, shp_path):
        """
        A shapefile in a read-only folder can be read but nothing can be written
        beside it, and several tools write their results next to the source. Copy
        the shapefile (all sidecar files) to a writable folder, once, and use that.
        """
        copies = self.__dict__.setdefault('_writable_copies', {})
        if shp_path in copies and os.path.exists(copies[shp_path]):
            return copies[shp_path]
        base = os.path.splitext(shp_path)[0]
        stem = self._safe_filename(os.path.basename(base))
        dest_base = os.path.join(self._fallback_dir(), stem)
        try:
            for ext in ('.shp', '.shx', '.dbf', '.prj', '.cpg', '.qpj', '.sbn', '.sbx'):
                if os.path.exists(base + ext):
                    shutil.copyfile(base + ext, dest_base + ext)
        except OSError as exc:
            self._warn(f"'{os.path.dirname(shp_path)}' is not writable and a working "
                       f"copy could not be made: {exc}", level=MSG_CRITICAL)
            return shp_path
        copies[shp_path] = dest_base + '.shp'
        self._warn(f"'{os.path.dirname(shp_path)}' is not writable, so a working copy "
                   f"was made in {os.path.dirname(dest_base)} and results will be saved there.",
                   level=MSG_INFO, duration=10)
        return copies[shp_path]

    def _make_output_path(self, qgis_layer, suffix):
        """
        A unique, safe .shp path for a result derived from `qgis_layer`:
        <folder>/<layer name><suffix>.shp, with a random number added if that
        file already exists.
        """
        base = os.path.join(self._output_directory(qgis_layer),
                            self._safe_filename(qgis_layer.name()) + self._safe_filename(suffix))
        path = base + ".shp"
        if os.path.exists(path):
            path = f"{base}_{random.randint(10000, 99999)}.shp"
        return path

    def _save_gdf(self, gdf, path):
        """
        Save a GeoDataFrame as a shapefile and return the path actually used.
        On failure, report the real reason and the path, retry in the temp folder,
        and only then raise (with the path in the message).
        """
        try:
            gdf.to_file(path, driver="ESRI Shapefile")
            return path
        except Exception as exc:
            self._warn(f"Could not save {path}: {exc}. Trying the temp folder instead.")
        alt = os.path.join(tempfile.gettempdir(), os.path.basename(path))
        try:
            gdf.to_file(alt, driver="ESRI Shapefile")
        except Exception as exc2:
            raise RuntimeError(f"could not save the result to {alt}: {exc2}") from exc2
        self._warn(f"Result saved to {alt} instead.", level=MSG_INFO)
        return alt

    def finalisePointSampler(self, gdf2, qgis_layer, name, param=0):
        layer_path = self._output_directory(qgis_layer)
        if(param==0):
            param2 = ''

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

        stem = self._safe_filename(qgis_layer.name())
        new_path = (
            layer_path
            + "/"
            + stem
            + "_"
            + self._safe_filename(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
                + "/"
                + stem
                + "_"
                + self._safe_filename(name)
                + param2
                + "_"
                + str(random_5_digit_integer)
                + ".shp"
            )
        # print("new_path_2:", new_path)
        # Write GeoPandas back to file
        if gdf2.empty:
            # An empty GeoDataFrame has no geometry to infer a type from and
            # gets written as a polyline layer; force an explicit Point schema.
            def _field_type(dtype):
                if pd.api.types.is_integer_dtype(dtype):
                    return "int"
                if pd.api.types.is_float_dtype(dtype):
                    return "float"
                return "str"
            schema = {
                "geometry": "Point",
                "properties": {
                    c: _field_type(gdf2[c].dtype)
                    for c in gdf2.columns if c != gdf2.geometry.name
                },
            }
            try:
                # pyogrio engine ignores `schema` but honours geometry_type
                gdf2.to_file(new_path, driver="ESRI Shapefile",
                             engine="pyogrio", geometry_type="Point")
            except Exception:
                gdf2.to_file(new_path, driver="ESRI Shapefile", schema=schema)
            self.iface.messageBar().pushMessage(
                "No points were retained - an empty point layer was created. "
                "Check that the polygon and point layers overlap.",
                level=MSG_WARNING,
                duration=15,
            )
        else:
            new_path = self._save_gdf(gdf2, new_path)

        # 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,
                                   count_fn=None):
        """
        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).

        count_fn(cell_size) -> number of output points, if it can be had without running
        the method (e.g. engine.count_cells); otherwise run_fn is run at every trial size.
        """
        area = (max_x - min_x) * (max_y - min_y)
        tol = 0.15

        def count(cs):
            if count_fn is not None:
                n = count_fn(cs)
            else:
                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 _restore_strike(self, gdf2, gdf_in, dip_dir_col):
        """
        The grid/Kent engines average dip direction. When the input field holds
        strike (Dip Direction box unticked) it was converted with
        dipdir = strike - 90, so convert the averaged values back, otherwise the
        output strikes come out rotated 90 degrees anticlockwise.
        """
        if (not self.dockwidget.checkBox_dip_dir.isChecked()
                and dip_dir_col in gdf2.columns):
            gdf2 = gdf2.copy()
            gdf2[dip_dir_col] = dipdir_to_strike(
                gdf2[dip_dir_col], gdf_in.get(dip_dir_col))
        return gdf2

    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

        # Optional kappa filter: blank (or 0) keeps every cell
        kappa_text = self.dockwidget.lineEdit_kent_kappa_min.text().strip()
        try:
            kappa_min = float(kappa_text) if kappa_text else None
        except ValueError:
            self._warn("Kent: Min kappa must be a number (or blank for no filter).")
            return
        if kappa_min is not None and kappa_min <= 0:
            kappa_min = None

        def kent_fn(*args, **kwargs):
            return engine.spherical_kent(*args, kappa_min=kappa_min, **kwargs)

        if self.dockwidget.checkBox_target_n.isChecked():
            target_n = self.dockwidget.spinBox_target_n.value()
            # Without a kappa filter every non-empty cell is output, so just count cells
            count_fn = None if kappa_min is not None else (
                lambda cs: engine.count_cells(gdf, min_x, max_x, min_y, max_y, cs))
            grid_size, _ = self._find_grid_size_for_target(
                kent_fn, gdf, min_x, max_x, min_y, max_y, target_n, count_fn=count_fn)
            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 = kent_fn(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 None and kappa_min is not None:
            self._warn(f"Kent: no cell reached Min kappa {kappa_min:g} - nothing written. "
                       f"Lower Min kappa or enlarge the grid size.")
        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})
            gdf2 = self._restore_strike(gdf2, gdf, dip_dir_col)
            name = ("grid_cell_kent" if kappa_min is None
                    else f"grid_cell_kent_kappa{kappa_min:g}")
            self.finalisePointSampler(gdf2, self.points_layer, name, 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

        # "Iterative": keep removing outliers (leave-one-out, repeated) instead of just one
        iterative = self.dockwidget.checkBox_kent_outlier_iter.isChecked()

        def outlier_fn(*args, **kwargs):
            return engine.outlier_removal(*args, iterative=iterative, **kwargs)

        if self.dockwidget.checkBox_target_n.isChecked():
            target_n = self.dockwidget.spinBox_target_n.value()
            # Every non-empty cell is output, whether or not outliers are removed from it
            grid_size, _ = self._find_grid_size_for_target(
                outlier_fn, gdf, min_x, max_x, min_y, max_y, target_n,
                count_fn=lambda cs: engine.count_cells(gdf, min_x, max_x, min_y, max_y, cs))
            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 = outlier_fn(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})
            gdf2 = self._restore_strike(gdf2, gdf, dip_dir_col)
            name = "grid_cell_kentOutlierIter" if iterative else "grid_cell_kentOutlier"
            self.finalisePointSampler(gdf2, self.points_layer, name, display_size)

    def carmichael(self):
        """Grid-cell Kent averaging with iterative multi-outlier removal (Carmichael 2016)."""
        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()
            # Carmichael outputs cells with more than 3 measurements
            grid_size, _ = self._find_grid_size_for_target(
                engine.outlier_carmichael, gdf, min_x, max_x, min_y, max_y, target_n,
                count_fn=lambda cs: engine.count_cells(gdf, min_x, max_x, min_y, max_y, cs,
                                                       min_points=4))
            display_size = grid_size * 110000 if is_geographic else grid_size
            self.dockwidget.lineEdit_grid_size_carm.setText(f"{display_size:.0f}")
        else:
            grid_size = float(self.dockwidget.lineEdit_grid_size_carm.text())
            display_size = grid_size
            if is_geographic:
                grid_size = grid_size / 110000

        tmpdir = tempfile.mkdtemp()
        try:
            file = engine.outlier_carmichael(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 None:
            self._warn("Carmichael: no cell had more than 3 measurements - nothing written. "
                       "Use a larger grid size.")
        else:
            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})
            gdf2 = self._restore_strike(gdf2, gdf, dip_dir_col)
            self.finalisePointSampler(gdf2, self.points_layer, "grid_cell_carmichael", 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=MSG_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=MSG_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_kent.isChecked():
            self.kent()
        elif self.dockwidget.radioButton_kent_outlier.isChecked():
            self.kentOutlier()
        elif self.dockwidget.radioButton_carmichael.isChecked():
            self.carmichael()
        elif self.dockwidget.radioButton_1o.isChecked():
            self.firstOrder()

    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:
                self._warn(f"T&P: fault 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:
                self._warn(f"T&P: fault cluster pre-compute failed: {exc}")
            finally:
                shutil.rmtree(tmpdir, ignore_errors=True)

        elif method == 'strat_offset':
            has_strat = ('max_off_lv' in gdf.columns or 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:
                    self._warn(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=MSG_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=MSG_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_kent.isChecked():
            point_method = 'spherical_kent'
        elif self.dockwidget.radioButton_1o.isChecked():
            point_method = 'firstorder'
        elif self.dockwidget.radioButton_carmichael.isChecked():
            point_method = 'outlier_carmichael'
        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._parse_codes(
            self.dockwidget.plainTextEdit_dyke_Codes.toPlainText())
        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_kent.text())
        except ValueError:
            grid_sz = 5000.0

        # First order needs the polygon boundaries as contacts, plus the
        # distance/angle limits; the T&P target count then picks the best-ranked
        # points within those limits.
        firstorder_kwargs = None
        if point_method == 'firstorder' and current_pt_gdf is not None:
            if current_polygon_gdf is None:
                self._warn("T&P with 1o Sampling needs a polygon layer for the contacts.")
                return
            try:
                fo_dist = float(self.dockwidget.lineEdit_1o_distance.text())
                fo_angle = float(self.dockwidget.lineEdit_1o_angle.text())
            except ValueError:
                self._warn("T&P with 1o Sampling: enter valid numbers for Distance and Angle.")
                return
            fo_contacts = current_polygon_gdf[['geometry']].copy()
            fo_contacts['geometry'] = current_polygon_gdf.boundary
            firstorder_kwargs = {
                'contact_gdf': fo_contacts,
                'dist_buffer': fo_dist,
                'angle_tol': fo_angle,
            }

        # 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

        # --- Report: one row per layer per step, written when the run finishes ---
        report_rows = []
        orig_counts = {
            'points': len(original_pt_gdf) if original_pt_gdf is not None else None,
            'faults': len(current_fault_gdf) if current_fault_gdf is not None else None,
            'polygons': len(current_polygon_gdf) if current_polygon_gdf is not None else None,
        }
        exponent = {'points': 1, 'faults': 2, 'polygons': 3}

        def add_report(step, ratio_val, step_ratio_val, layer, method, n_in, n_out, info):
            on0 = orig_counts[layer]
            report_rows.append({
                'step': step,
                'OS/TS ratio': ratio_val,
                'step ratio': step_ratio_val,
                'layer': layer,
                'method': method,
                'input count': n_in,
                'target count (from original)': topfer_count(on0, ratio_val, exponent[layer]),
                'output count': n_out,
                'parameters used': '; '.join(
                    f"{k}={v:.6g}" if isinstance(v, float) else f"{k}={v}"
                    for k, v in info.items() if k != 'method'),
            })

        # --- 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")
            add_report('baseline', 1.0, 1.0, 'points', 'none (all kept)',
                       len(original_pt_gdf), len(original_pt_gdf), {})

        if current_fault_gdf is not None:
            ln_full = scale_lines_tp(current_fault_gdf, 1.0, fault_method)
            add_report('baseline', 1.0, 1.0, 'faults', fault_method,
                       len(current_fault_gdf), len(ln_full), {})
            out_path = self._save_gdf(ln_full, self._make_output_path(fault_layer, "_tp_1"))
            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:
            out_path = self._save_gdf(current_polygon_gdf,
                                      self._make_output_path(polygon_layer, "_tp_1"))
            poly_layer = QgsVectorLayer(out_path, f"{polygon_layer.name()}_tp_1", "ogr")
            self._add_layer_with_style(poly_layer, polygon_layer)
            add_report('baseline', 1.0, 1.0, 'polygons', 'none (all kept)',
                       len(current_polygon_gdf), len(current_polygon_gdf), {})

        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_info = {}
                pt_scaled = scale_points_tp(
                    original_pt_gdf, iter_ratio, pt_engine, point_method,
                    pt_dip_col, pt_dip_dir_col, grid_sz,
                    firstorder_kwargs=firstorder_kwargs, info=pt_info,
                    strike_input=not self.dockwidget.checkBox_dip_dir.isChecked())
                add_report(f"ratio {ratio_tag}", iter_ratio, iter_ratio, 'points',
                           point_method, len(original_pt_gdf), len(pt_scaled), pt_info)
                tn_expected = topfer_count(len(original_pt_gdf), iter_ratio, x=1)
                if point_method == 'firstorder' and len(pt_scaled) < tn_expected:
                    self._warn(
                        f"T&P ratio {iter_ratio}: only {len(pt_scaled)} points pass the 1o "
                        f"distance/angle limits (target {tn_expected}). Loosen the limits for more.",
                        duration=10)

                # 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_info = {}
                n_ln_in = len(current_fault_gdf)
                ln_scaled = scale_lines_tp(current_fault_gdf, step_ratio, fault_method,
                                           info=ln_info)
                add_report(f"ratio {ratio_tag}", iter_ratio, step_ratio, 'faults',
                           fault_method, n_ln_in, len(ln_scaled), ln_info)
                current_fault_gdf = ln_scaled
                out_path = self._save_gdf(
                    ln_scaled, self._make_output_path(fault_layer, f"_tp_{ratio_tag}"))
                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:
                poly_info = {}
                n_poly_in = len(current_polygon_gdf)
                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),
                        info=poly_info,
                    )
                    current_polygon_gdf = poly_scaled
                except Exception as exc:
                    self._warn(f"T&P polygon scaling failed at ratio {iter_ratio}; "
                               f"unscaled polygons written: {exc}")
                    poly_info['note'] = f"FAILED, unscaled polygons written: {exc}"
                    poly_scaled = current_polygon_gdf.copy()
                add_report(f"ratio {ratio_tag}", iter_ratio, step_ratio, 'polygons',
                           'merge smallest into best-matching neighbour',
                           n_poly_in, len(poly_scaled), poly_info)

                out_path = self._save_gdf(
                    poly_scaled, self._make_output_path(polygon_layer, f"_tp_{ratio_tag}"))
                poly_layer = QgsVectorLayer(
                    out_path,
                    f"{polygon_layer.name()}_tp_{ratio_tag}", "ogr")
                self._add_layer_with_style(poly_layer, polygon_layer)

        # --- Write the T&P report (CSV) next to the first available input layer ---
        settings = {
            'target OS/TS ratio': ratio,
            'increment': increment,
            'points layer': self.points_layer.name() if original_pt_gdf is not None else 'none',
            'point method': point_method if original_pt_gdf is not None else 'n/a',
            'fault layer': fault_layer.name() if current_fault_gdf is not None else 'none',
            'fault method': fault_method if current_fault_gdf is not None else 'n/a',
            'polygon layer': polygon_layer.name() if polygon_layer else 'none',
            'polygon priority fields': ', '.join(
                c for c in [strat1, strat2, strat3, strat4, lithoname] if c) or 'none',
            'polygon node tolerance': dist_thresh,
        }
        if point_method == 'firstorder' and firstorder_kwargs:
            settings['1o distance limit'] = firstorder_kwargs['dist_buffer']
            settings['1o angle limit'] = firstorder_kwargs['angle_tol']
        elif point_method in ('gridcell_average', 'spherical_kent', 'outlier_removal',
                              'outlier_carmichael'):
            settings['grid size starting guess'] = grid_sz
        report_rows.insert(0, {
            'step': 'settings', 'OS/TS ratio': '', 'step ratio': '', 'layer': 'all',
            'method': '', 'input count': '', 'target count (from original)': '',
            'output count': '',
            'parameters used': '; '.join(f"{k}={v}" for k, v in settings.items()),
        })
        try:
            ref_layer = next(l for l in (self.points_layer if original_pt_gdf is not None else None,
                                         fault_layer if current_fault_gdf is not None else None,
                                         polygon_layer if current_polygon_gdf is not None else None)
                             if l is not None)
            report_dir = os.path.dirname(self._resolved_source(ref_layer))
            report_path = os.path.join(
                report_dir, f"TP_scaling_report_{time.strftime('%Y%m%d_%H%M%S')}.csv")
            with open(report_path, 'w', newline='', encoding='utf-8-sig') as fh:
                writer = csv.DictWriter(fh, fieldnames=list(report_rows[0].keys()))
                writer.writeheader()
                writer.writerows(report_rows)
            self.iface.messageBar().pushMessage(
                f"T&P scaling complete. Report written to {report_path}",
                level=MSG_SUCCESS, duration=15)
            print(f"T&P report: {report_path}")
        except Exception as exc:
            self._warn(f"T&P scaling finished but the report could not be written: {exc}")

    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:
                    self._warn("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)
                problem = dip_field_problem(orientation_gdf[dip_col], dip_col)
                if problem:
                    self._warn(problem + " Nothing was processed.", level=MSG_CRITICAL, duration=20)
                    return

                # The engine expects dip direction. If the field holds strike,
                # convert it exactly as setUpPointSampler does (the engine's
                # strike = dipdir - 90 then recovers the original strike mod 180).
                eff_dipdir_col = dip_dir_col
                if not self.dockwidget.checkBox_dip_dir.isChecked():
                    orientation_gdf['_eff_dipdir'] = (
                        orientation_gdf[dip_dir_col].astype(float) - 90) % 360
                    eff_dipdir_col = '_eff_dipdir'

                # Create an instance of SubsamplingEngine
                subsampler = SubsamplingEngine(dip=dip_col, dipdir=eff_dipdir_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)
                if output_gdf is not None and '_eff_dipdir' in output_gdf.columns:
                    output_gdf = output_gdf.drop(columns=['_eff_dipdir'])

                if output_gdf is None or output_gdf.empty:
                    print("First Order: no measurements retained - check the polygon "
                          "and point layers overlap and the distance/angle settings.")
                    # the message-bar warning is raised by finalisePointSampler
                    if output_gdf is None:
                        return

                # 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._parse_codes(
                    self.dockwidget.plainTextEdit_dyke_Codes.toPlainText())
                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

                    new_path = self._make_output_path(
                        self.polygon_layer, "_min_area_" + parameter)

                    # 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 and not dyke_codes:
                        self._warn("Ignore field is set but the Ignore Codes box is empty - "
                                   "no polygons will be triangulated.")
                    if dyke_index != 0 and dyke_codes:
                        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)
                        self._report_ignore_codes(triangulator, dyke_codes, gdf)

                    # 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,
                    )
                    new_path = self._save_gdf(output_gdf, new_path)

                    # 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=MSG_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°)
            )

            new_path = self._make_output_path(self.polyline_layer, "_merged")

            gdf = gpd.read_file(self._resolved_source(self.polyline_layer))
            output_gdf = merger.process_shapefile(gdf)
            new_path = self._save_gdf(output_gdf, new_path)

            # 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=MSG_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:
                self._warn("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=MSG_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"
            )

            safe_layer_name = self._safe_filename(self.polyline_layer.name())
            best_n = filter_obj.example_manual_clusters(
                outPath,
                safe_layer_name,
                azimuth_field="endpt_az",
            )
            if best_n is None:
                self.iface.messageBar().pushMessage(
                    "Could not cluster faults - too few faults with distinct orientations (see console).",
                    level=MSG_WARNING,
                    duration=15,
                )
                return
            # the clustering step writes beside outPath, using the same safe name
            new_path = os.path.join(
                os.path.dirname(outPath), safe_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:
            self._warn("Map simplification: select both a polygon layer and a fault layer.")
            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_kent.setToolTip(
            "Grid Cell Averaging of points using Kent statistics.\n"
            "Output points sit at the mean position of each cell's measurements, with the vector-mean\n"
            "orientation plus count, Kent kappa (concentration) and beta (ovalness) attributes."
        )
        self.dockwidget.lineEdit_grid_size_kent.setToolTip(
            "Grid size in layer units for cell averaging using Kent statistics"
        )
        self.dockwidget.lineEdit_kent_kappa_min.setToolTip(
            "Kent only. Keep a cell only if its Kent kappa (concentration) is at least this value.\n"
            "Drops cells with a single measurement and cells whose orientations are scattered.\n"
            "Blank (or 0) keeps every cell. The thesis models used 5."
        )
        self.dockwidget.radioButton_carmichael.setToolTip(
            "Grid Cell Averaging using Kent statistics, removing several outliers per cell\n"
            "(Carmichael & Ailleres 2016). Points are ranked by angle from the cell mean and\n"
            "removed one at a time while Kent kappa is recomputed; everything up to the\n"
            "biggest jump in kappa is treated as outliers. Cells with 3 or fewer points are skipped."
        )
        self.dockwidget.checkBox_kent_outlier_iter.setToolTip(
            "Kent Outlier only. Remove several outliers per cell instead of one: points are taken\n"
            "out in order of how much their removal raises Kent kappa (up to half the cell), and the\n"
            "sequence is cut at its biggest kappa jump - only if that jump stands out (5x the\n"
            "typical step), otherwise nothing is removed. Differs from Carmichael, which orders\n"
            "points by angle from the cell mean. Not applied by the Töpfer & Pillewizer tool."
        )
        self.dockwidget.lineEdit_grid_size_carm.setToolTip(
            "Grid size in layer units for Carmichael outlier averaging"
        )
        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"
        )

        target_n_tip = (
            "Kent, Kent Outlier and Carmichael only.\n"
            "Instead of using the Grid Size box, search for the grid cell size that gives "
            "approximately N output points (within ±15%). The grid size found is written back "
            "into the Grid Size box of the selected method.\n"
            "Several trial runs are needed, so this is slower. Ignored by Stochastic and 1o Sampling."
        )
        self.dockwidget.checkBox_target_n.setToolTip(target_n_tip)
        self.dockwidget.spinBox_target_n.setToolTip(
            "Desired number of output points (used when 'Auto-fit to N pts' is ticked). "
            "The result is approximate (±15%) and cannot exceed the number of occupied grid cells."
        )
        self.dockwidget.radioButton_1o.setToolTip(
            "First order retention of points based on distance and angle to contacts.\n"
            "Requires a polygon layer (selected in the map polygons box) that overlaps the points - "
            "its boundaries are used as the 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 ignore codes")
        self.dockwidget.plainTextEdit_dyke_Codes.setToolTip(
            "Comma separated list of ignore 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(
                POINT_LAYER
            )

            self.dockwidget.mMapLayerComboBox_maps_polygons.setFilters(
                POLYGON_LAYER
            )

            self.dockwidget.mMapLayerComboBox_fault_polylines.setFilters(
                LINE_LAYER
            )
            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._guard(self.subsamplePoints, "Point subsampling"))

            self.dockwidget.pushButton_minPolyArea.clicked.connect(
                self._guard(self.minPolyArea, "Polygon cleaning"))

            self.dockwidget.pushButton_merge_segments.clicked.connect(
                self._guard(self.mergeSegments, "Fault segment merging")
            )
            self.dockwidget.pushButton_process_faults.clicked.connect(
                self._guard(self.processFaults, "Fault processing"))
            self.dockwidget.pushButton_tp_run.clicked.connect(
                self._guard(self.topferPillewizer, "Töpfer & Pillewizer scaling"))
            self.dockwidget.pushButton_simplifyMap.clicked.connect(
                self._guard(self.simplifyMap, "Map simplification")
            )
            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)