# -*- coding: utf-8 -*-
"""
/***************************************************************************
 SGTool
                                 A QGIS plugin
 Simple Potential Field Processing
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2024-11-17
        git sha              : $Format:%H$
        copyright            : (C) 2024 by Mark Jessell
        email                : mark.jessell@uwa.edu.au
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.QtGui import QIcon, QDesktopServices, QValidator
from qgis.PyQt.QtWidgets import QAction, QDockWidget, QFileDialog, QMessageBox
from qgis.core import (
    Qgis,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsVectorLayer,
    QgsProject,
    QgsRasterLayer,
    QgsSingleBandGrayRenderer,
    QgsSingleBandPseudoColorRenderer,
    QgsColorRampShader,
    QgsRasterShader,
    QgsStyle,
    QgsCptCityColorRamp,
    QgsFeature,
    QgsField,
    QgsProcessingFeedback,
    QgsProcessingFeatureSourceDefinition,
    QgsFeatureRequest,
    QgsWkbTypes,
    QgsProject,
    QgsGeometry,
    QgsFields,
    QgsPointXY,
    QgsMapLayerProxyModel,
    QgsApplication,
    QgsVectorFileWriter,
    QgsVectorLayer,
    QgsCoordinateTransformContext,
    QgsMapLayerType,
)
from qgis.gui import QgsMapToolEmitPoint

from qgis.PyQt.QtCore import (
    QSettings,
    QTranslator,
    QCoreApplication,
    QFileInfo,
    QVariant,
    Qt,
    QUrl,
    QLocale,
)


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

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

    # QMessageBox buttons
    QMessageBox_Ok = QMessageBox.StandardButton.Ok
    QMessageBox_Cancel = QMessageBox.StandardButton.Cancel
    QMessageBox_Yes = QMessageBox.StandardButton.Yes
    QMessageBox_No = QMessageBox.StandardButton.No

except AttributeError:
    # Qt5 detected
    QT6 = False

    # Qt5 style enums
    RightDockWidgetArea = Qt.RightDockWidgetArea
    LeftDockWidgetArea = Qt.LeftDockWidgetArea
    TopDockWidgetArea = Qt.TopDockWidgetArea
    BottomDockWidgetArea = Qt.BottomDockWidgetArea

    # QMessageBox buttons
    QMessageBox_Ok = QMessageBox.Ok
    QMessageBox_Cancel = QMessageBox.Cancel
    QMessageBox_Yes = QMessageBox.Yes
    QMessageBox_No = QMessageBox.No

# from PyQt5.QtGui import QValidator

import re
import os.path
import numpy as np
from scipy.spatial import cKDTree
from scipy import interpolate
from scipy.interpolate import interp1d
import tempfile
from datetime import datetime
from pyproj import Transformer
import processing
from osgeo import gdal, osr
import platform
import importlib
import subprocess
import os

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

# Import the code for the DockWidget
from .SGTool_dockwidget import SGToolDockWidget
from .calcs.GeophysicalProcessor import GeophysicalProcessor
from .calcs.Cooper_Cowan_rtpvariable import rtpvariable as cooper_rtpvariable
from .calcs.geosoft_grid_parser import *
from .calcs.PSplot import PowerSpectrumDock
from .calcs.ConvolutionFilter import ConvolutionFilter

# from .calcs.ConvolutionFilter import OddPositiveIntegerValidator
from .calcs.GridData_no_pandas import QGISGridData
from .calcs.SpatialStats import SpatialStats
from .calcs.WTMM import WTMM
from .calcs.PCAICA import PCAICA
from .calcs.SG_Util import SG_Util
from .calcs.igrf.igrf_utils import igrf_utils as IGRF
from .calcs.aseggdf2parser import AsegGdf2Parser

# from .calcs.euler.euler_python_optimised import euler_deconv_opt
from .calcs.euler.euler_python import euler_deconv
from .calcs.euler.estimates_statistics import window_stats
from .calcs.mrvbf import mrvbf as calc_mrvbf
from .calcs.saga_mba_gridding import mba_gridding


class RGBPickerMapTool(QgsMapToolEmitPoint):
    """One-shot map tool: emits a point on click then deactivates itself."""

    def __init__(self, canvas):
        super().__init__(canvas)
        try:
            self.setCursor(Qt.CursorShape.CrossCursor)
        except AttributeError:
            self.setCursor(Qt.CrossCursor)

    def canvasPressEvent(self, event):
        pass  # consume press; do not propagate

    def canvasReleaseEvent(self, event):
        # Emit signal ourselves, do NOT call super() — prevents QGIS propagating
        # this click to the identify handler or any other default behaviour.
        point = self.toMapCoordinates(event.pos())
        self.canvasClicked.emit(point, event.button())


class SGTool:
    """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 — try full locale code first (e.g. pt_BR), then 2-char fallback (e.g. pt)
        full_locale = str(QSettings().value("locale/userLocale"))
        locale_path = os.path.join(
            self.plugin_dir, "i18n", "SGTool_{}.qm".format(full_locale)
        )
        if not os.path.exists(locale_path):
            locale_path = os.path.join(
                self.plugin_dir, "i18n", "SGTool_{}.qm".format(full_locale[0:2])
            )

        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("&SGTool")
        # TODO: We are going to let the user set this up in a future iteration
        self.toolbar = self.iface.addToolBar("SGTool")
        self.toolbar.setObjectName("SGTool")

        # print "** INITIALIZING SGTool"

        self.pluginIsActive = False
        self.dlg = None
        self.last_directory = 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("SGTool", 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/SGTool/icon.png"
        self.add_action(
            icon_path,
            text=self.tr("SGTool"),
            callback=self.run,
            parent=self.iface.mainWindow(),
        )

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

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

        # print "** CLOSING SGTool"

        # disconnects
        self.dlg.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.dlg = None

        self.pluginIsActive = False

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

        # print "** UNLOAD SGTool"

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

    def define_tips(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.
        Tooltips are set for:
        - Map layer combo boxes
        - Buttons for file selection and processing
        - Checkboxes for various processing options
        - Input fields for parameters like azimuth, inclination, declination, etc.
        - Dropdowns for selecting options like derivative direction, continuation type, etc.
        - Spin boxes for numerical inputs like cell size, window size, etc.
        - Group boxes for organizing related functionalities
        - Labels providing additional context for specific operations
        The tooltips cover a wide range of geophysical processing tasks,
        including filtering, reduction to pole/equator, derivative calculations,
        analytic signal, continuation, bandpass filtering, spatial statistics,
        and more. They also include explanations for advanced operations like
        wavelet transforms, DTM classification, and geotiff normalization.
        This method enhances the user experience by providing inline documentation
        for the plugin's features, making it easier for users to navigate and
        utilize the available tools effectively.
        """

        self.dlg.mMapLayerComboBox_selectGrid.setToolTip(
            self.tr("Layer selected for processing")
        )
        self.dlg.mMapLayerComboBox_selectGrid_Conv.setToolTip(
            self.tr("Layer selected for grid processing")
        )
        self.dlg.pushButton_2_selectGrid.setToolTip(
            self.tr("Load new file for processing\n*.grd *.tif *.ers *.grv *.mag, for other formats use normal QGIS grid loader")
        )
        self.dlg.checkBox_3_DirClean.setToolTip(
            self.tr("Filter (DirCos + Butterworth) remove a specific direction and wavelength,\nUseful for filtering flight line noise")
        )

        self.dlg.lineEdit_3_azimuth.setToolTip(
            self.tr("Azimuth of high frequency noise to be filtered (degrees clockwise from North)")
        )

        self.dlg.checkBox_4_RTE_P.setToolTip(
            self.tr("Reduction to pole or equator\nThe reduction to the pole (RTP), differential RTP or to Equator (RTE) is a process in geophysics\nwhere magnetic data are transformed to look as though\n they were measured at the magnetic pole/equator\nCorrects the asymmetry of magnetic anomalies caused by\n the Earth's field, making them appear directly above their sources\n Differential RTP (code derived from Cooper & Cowan 2005) calculates RTP with spatially varying inclination and declination")
        )
        self.dlg.pushButton_4_calcIGRF.setToolTip(
            self.tr("Calculate IGRF Inclination & Declination based on centroid of selected grid and specified survey height and date\nIf a grid originated from a Noddy calculation, the info is read from the geotiff metadata")
        )
        self.dlg.comboBox_3_rte_p_list.setToolTip(
            self.tr("RTP: Reduction to Pole (high mag latitudes >20°)\nDiff. RTP: Variable (Differential) RTP using Cooper/Cowan Taylor series — accounts for spatial variation in inc/dec across the survey area\nRTE: Reduction to Equator (low mag latitudes <20°)")
        )
        self.dlg.lineEdit_6_inc.setToolTip(
            self.tr("Manually define magnetic inclination [degrees from horizontal]")
        )
        self.dlg.lineEdit_5_dec.setToolTip(
            self.tr("Manually define magnetic declination [degrees clockwise from North]")
        )
        self.dlg.lineEdit_6_int.setToolTip(self.tr("Survey intensity in nT"))
        self.dlg.dateEdit.setToolTip(self.tr("Survey date (1900-2030)"))
        self.dlg.checkBox_4_PGrav.setToolTip(
            self.tr("Vertical Integration:\nWhen applied to RTE/P result converts magnetic anomalies into gravity-like anomalies (i.e. same decay with distance from source) for comparison or joint interpretation\nAlso good for stitched grids with very different line spacing.\nRequires a metre-based projection")
        )

        self.dlg.checkBox_5_regional.setToolTip(
            self.tr("Remove regional (RR) based on 1st or 2nd order polynomial ")
        )

        self.dlg.checkBox_6_derivative.setToolTip(
            self.tr("Calculate derivate (d+power+direction) parallel to x, y or z\nHighlights near-surface/short-wavelength features")
        )

        self.dlg.comboBox_derivDirection.setToolTip(self.tr("Select derivative direction"))
        self.dlg.lineEdit_9_derivePower.setToolTip(self.tr("Power of derivative"))
        self.dlg.checkBox_7_tiltDerivative.setToolTip(
            self.tr("Tilt Derivative (TD)\nIt is often applied to magnetic or gravity data to enhance edges and detect shallow sources\nTends to overconnect structural features")
        )

        self.dlg.checkBox_8_analyticSignal.setToolTip(
            self.tr("Analytic Signal (AS)\nIt combines horizontal and vertical derivatives to highlight anomaly edges and amplitude variations, independent of direction")
        )

        self.dlg.checkBox_9_continuation.setToolTip(
            self.tr("Upward or downward continuation\nUpward Continuation (UC) data by continuing it to a higher altitude, attenuating high-frequency noise and shallow features\nDownward Continuation (DC) enhances shallow or high-frequency anomalies by continuing the field to a lower altitude")
        )

        self.dlg.comboBox_2_continuationDirection.setToolTip(
            self.tr("Select direction of continuation")
        )
        self.dlg.lineEdit_10_continuationHeight.setToolTip(
            self.tr("Select amount of continuation [m only]")
        )
        self.dlg.checkBox_10_bandPass.setToolTip(
            self.tr("Band pass filter (BP)\nIsolates specific wavelength features.")
        )

        self.dlg.lineEdit_12_bandPassLow.setToolTip(
            self.tr("Low wavelength cutoff [m or other length unit]")
        )
        self.dlg.lineEdit_11_bandPassHigh.setToolTip(
            self.tr("High wavelength cutoff [m or other length unit]")
        )
        self.dlg.lineEdit_3_HLP_width.setToolTip(
            self.tr("Width of cosine rolloff [m or other length unit]\nStart with cutoff value and increase to reduce ringing\n0 is step cutoff")
        )
        self.dlg.lineEdit_3_BP_width.setToolTip(
            self.tr("Width of cosine rolloff [m or other length unit]\nStart with cutoff value and increase to reduce ringing\n0 is step cutoff")
        )

        self.dlg.checkBox_10_freqCut.setToolTip(
            self.tr("High or Low pass filter\nIsolates specific short wavelength (HP) or long wavelength (LP) features.")
        )

        self.dlg.comboBox_2_FreqCutType.setToolTip(self.tr("Cut off type"))
        self.dlg.lineEdit_12_FreqPass.setToolTip(
            self.tr("Cutoff wavelength [m or other length unit]")
        )
        self.dlg.checkBox_11_1vd_agc.setToolTip(
            self.tr("Automatic Gain Control (AGC) or Amplitude Normalisation\nHighlights short wavelength/low amplitude features")
        )

        self.dlg.lineEdit_13_agc_window.setToolTip(self.tr("Window size for normalisation"))
        self.dlg.pushButton_3_applyProcessing.setToolTip(
            self.tr("Apply selected processing steps in parallel to selected grid")
        )
        self.dlg.pushButton_3_applyProcessing_Conv.setToolTip(
            self.tr("Apply selected processing steps in parallel to selected grid")
        )
        self.dlg.lineEdit_13_max_buffer.setToolTip(
            self.tr("Maximum buffer to be applied to grid to reduce edge effects")
        )
        self.dlg.checkBox_11_tot_hz_grad.setToolTip(
            self.tr("Total Horizontal Gradient Calculation (THG)")
        )
        self.dlg.pushButton_rad_power_spectrum.setToolTip(
            self.tr("Provides pop-up display of grid plus Radial Averaged Power Spectrum (needs testing!)")
        )

        self.dlg.checkBox_Mean.setToolTip(
            self.tr("Mean of values around central pixel\nSmooths data")
        )

        self.dlg.checkBox_Median.setToolTip(
            self.tr("Median of values around central pixel\nRemoves high frequency noise")
        )

        self.dlg.checkBox_Gaussian.setToolTip(self.tr("Gaussian smoothing of image"))

        self.dlg.checkBox_Directional.setToolTip(
            self.tr("Directional enhancement\nHighlights high frequency data in a particular direction")
        )

        self.dlg.pushButton_selectPoints.setToolTip(
            self.tr("Select CSV, DAT (ASEG-GDF2 or legacy DAT) or XYZ format points file\nLegacy DAT format is quite variable so columns may be offset from field name!")
        )
        self.dlg.comboBox_grid_x.setToolTip(
            self.tr("Define X coordinate column (for csv & dat files)")
        )
        self.dlg.comboBox_grid_y.setToolTip(
            self.tr("Define Y coordinate column (for csv & dat files)")
        )
        self.dlg.mQgsProjectionSelectionWidget.setToolTip(
            self.tr("DEfine Coordinate System of point data")
        )
        self.dlg.checkBox_load_tie_lines.setToolTip(
            self.tr("For xyz format files only, optionally load tie lines")
        )
        self.dlg.pushButton_load_point_data.setToolTip(
            self.tr("Load points file and convert to layer\nWith polyline layer of lines for xyz format files")
        )
        self.dlg.mMapLayerComboBox_selectGrid_3.setToolTip(
            self.tr("Select from currently loaded points layers for gridding")
        )
        self.dlg.comboBox_select_grid_data_field.setToolTip(self.tr("Select field to grid"))
        self.dlg.doubleSpinBox_cellsize.setToolTip(self.tr("Define cell size in layer units"))
        self.dlg.pushButton_idw_2.setToolTip(
            self.tr("Perform Inverse Distance Weighting (IDW) gridding\nOpens Standard QGIS Dialog")
        )
        self.dlg.pushButton_bspline_3.setToolTip(
            self.tr("Perform Multilevel B-Spline (MBA) gridding of the selected point field\n"
                    "into a temporary grid (local NumPy translation of SAGA's grid_spline tool)")
        )
        self.dlg.lineEdit_bspline_epsilon.setToolTip(
            self.tr("MBA threshold error (EPSILON): stop refining once every point residual is below this value")
        )
        self.dlg.spinBox_bspline_levels.setToolTip(
            self.tr("MBA maximum number of refinement levels (LEVEL_MAX)")
        )
        self.dlg.checkBox_bspline_ignore.setToolTip(
            self.tr("Drop points whose data value is below the threshold before gridding\n"
                    "(useful for removing nodata sentinels such as -99999 or -999999.9)")
        )
        self.dlg.lineEdit_bspline_ignore.setToolTip(
            self.tr("Points with a value less than this are excluded from B-Spline gridding")
        )
        self.dlg.label_51.setToolTip(
            self.tr("Number of cells in x & y directions based on spatial extent of points and Cell Size")
        )
        self.dlg.pushButton_2_selectGrid_RGB.setToolTip(
            self.tr("Select RGB image that you want to convert to a monotonic grayscale image")
        )
        self.dlg.textEdit_2_colour_list.setToolTip(
            self.tr("Comma separated list of CSS colours\nOR a set of comma seperated RBG triplets")
        )
        self.dlg.groupBox_7.setToolTip(
            self.tr("1) Load a RGB raster image,\n2) Define a Look Up Table by defining a comma separated sequence of colours using CSS colour names OR a set of comma seperated RBG triplets and\n3) Convert to monotonically increasing greyscale image\n\nDo not use if any shading has been applied to the image!")
        )
        self.dlg.mQgsDoubleSpinBox_LUT_min.setToolTip(
            self.tr("Define min and max values for rescaling of grid values")
        )
        self.dlg.mQgsDoubleSpinBox_LUT_max.setToolTip(
            self.tr("Define min and max values for rescaling of grid values")
        )
        self.dlg.pushButton_CSSS_Colours.setToolTip(
            self.tr("See full suite of CSS Colours (requires network connection)")
        )
        self.dlg.lineEdit.setToolTip(self.tr("Example colour sequence, can be copy pasted"))

        self.dlg.spinBox_levels.setToolTip(self.tr("Number of levels"))
        self.dlg.doubleSpinBox_base.setToolTip(self.tr("Lowest height to worm"))
        self.dlg.doubleSpinBox_inc.setToolTip(self.tr("Increment in metres between levels "))
        self.dlg.groupBox_8.setToolTip(self.tr("Create csv file of worms using bsdwormer code"))
        self.dlg.mMapLayerComboBox_selectGrid_worms.setToolTip(
            self.tr("Grid to be analysed by wavelet transform")
        )

        self.dlg.checkBox_NaN.setToolTip(self.tr("Threshold background values to NaN"))
        self.dlg.radioButton_NaN_Above.setToolTip(
            self.tr("Threshold background values above set value to NaN")
        )
        self.dlg.radioButton_NaN_Below.setToolTip(
            self.tr("Threshold background values below set value to NaN")
        )
        self.dlg.radioButton_NaN_Both.setToolTip(
            self.tr("Threshold background values between set values to NaN")
        )
        self.dlg.doubleSpinBox_NaN_Above.setToolTip(self.tr("Upper threshold value"))

        self.dlg.checkBox_worms_shp.setToolTip(
            self.tr("Convert worms to polyline shapefile\n(Can be slow and best to start worms at >=2000m)")
        )

        self.dlg.pushButton_select_normalise_in.setToolTip(
            self.tr("Directory with geotiffs to be normalised")
        )

        self.dlg.pushButton_select_normalise_out.setToolTip(
            self.tr("Directory where normalised geotiffs will be stored")
        )

        self.dlg.pushButton_normalise.setToolTip(
            self.tr("Normalise geotiffs:\n1. First order regional removed\n2. Normalise standard deviation using alpahbetical first file in input directory as reference\n3. Remove mean")
        )
        self.dlg.groupBox_10.setToolTip(
            self.tr("Normalise directory of geotiffs\nOnce normalised the QGIS merge tool produced a reasonable stitch of the grids\nAssumes same processing level for grids\nAssumes flight heights have been normalised by continuation\nMerge uses alpahbetical first file to define cell size\nAll grids in merge have to be same projection")
        )

        self.dlg.radioButton_normalise_1st.setToolTip(
            self.tr("1st order (flat plane) regional removed from grid")
        )

        self.dlg.radioButton_normalise_2nd.setToolTip(
            self.tr("2nd order polynomial regional removed from grid (slower)")
        )

        self.dlg.checkBox_polygons.setToolTip(self.tr("Create grid outline polygon(s) layer"))

        self.dlg.label_26.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )
        self.dlg.label_28.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )
        self.dlg.label_29.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )
        self.dlg.label_32.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )
        self.dlg.label_30.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )
        self.dlg.label_25.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )
        self.dlg.label_40.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )
        self.dlg.label_41.setToolTip(
            self.tr("Functions preceded by dot points should be calculated on RTE or RTP mag data")
        )

        self.dlg.checkBox_SS_Min.setToolTip(
            self.tr("Calculate minimum of values around central pixel")
        )
        self.dlg.checkBox_SS_Max.setToolTip(
            self.tr("Calculate maximum of values around central pixel")
        )
        self.dlg.checkBox_SS_StdDev.setToolTip(
            self.tr("Calculate standard deviation of values around central pixel\nWill be slow for larger grids and window sizes!!")
        )
        self.dlg.checkBox_SS_Variance.setToolTip(
            self.tr("Calculate variance of values around central pixel\nWill be slow for larger grids and window sizes!!")
        )
        self.dlg.checkBox_SS_Skewness.setToolTip(
            self.tr("Calculate skewness of values around central pixel\nWill be VERY slow for larger grids and window sizes!!")
        )
        self.dlg.checkBox_SS_Kurtosis.setToolTip(
            self.tr("Calculate kurtosis of values around central pixel\nWill be VERY slow for larger grids and window sizes!!")
        )
        self.dlg.checkBox_SS_Anisotropy.setToolTip(
            self.tr("Calculate local anisotropy magnitude (0–1) and dominant orientation (degrees)\nReturns two layers: _SS_AnisoMag and _SS_AnisoOrient")
        )
        self.dlg.checkBox_SS_ChainLength.setToolTip(
            self.tr("Score each pixel by the length of the oriented chain it belongs to\n(discrete 8-connected, fast)")
        )
        self.dlg.checkBox_SS_Streamline.setToolTip(
            self.tr("Score each pixel by sub-pixel forward+backward path length along the orientation field\n(handles curves; slower than Chain Length)")
        )
        self.dlg.doubleSpinBox_SS_AnisoThresh.setToolTip(
            self.tr("Minimum anisotropy magnitude for a pixel to be considered part of a feature (0–1)")
        )
        self.dlg.doubleSpinBox_SS_AngleTol.setToolTip(
            self.tr("Maximum orientation difference (degrees) between adjacent pixels to be chained")
        )
        self.dlg.spinBox_SS_MaxSteps.setToolTip(
            self.tr("Maximum number of sub-pixel steps per direction for Streamline Length")
        )
        self.dlg.spinBox_SS_SearchRadius.setToolTip(
            self.tr("Chain Length: how many pixels ahead to cast the orientation ray when searching for the next connected pixel (1=8-neighbours only, 2–3 recommended for oblique features)")
        )

        self.dlg.lineEdit_SS_Window.setToolTip(
            self.tr("Size of window for calculation of spatial statistics")
        )

        self.dlg.checkBox_MRVBF.setToolTip(
            self.tr("Multiresolution Valley Bottom Flatness (MRVBF) and Ridge Top\n"
                    "Flatness (MRRTF) after Gallant & Dowling 2003. Writes _MRVBF,\n"
                    "_MRRTF and _MRVBF_slope grids plus a 3-band (MRRTF, MRVBF,\n"
                    "Slope) composite. Requires a projected raster in metres.")
        )
        self.dlg.lineEdit_MRVBF_tslope.setToolTip(
            self.tr("Initial threshold for slope (percent) at the finest resolution")
        )
        self.dlg.lineEdit_MRVBF_tpctlv.setToolTip(
            self.tr("Threshold for lowness percentile (valley bottom flatness)")
        )
        self.dlg.lineEdit_MRVBF_tpctlr.setToolTip(
            self.tr("Threshold for upness percentile (ridge top flatness)")
        )
        self.dlg.lineEdit_MRVBF_pslope.setToolTip(
            self.tr("Shape parameter for the slope transformation")
        )
        self.dlg.lineEdit_MRVBF_ppctl.setToolTip(
            self.tr("Shape parameter for the elevation percentile transformation")
        )
        self.dlg.lineEdit_MRVBF_maxres.setToolTip(
            self.tr("Maximum resolution as a percentage of the grid diagonal")
        )

        self.dlg.lineEdit_3_DC_wavelength.setToolTip(
            self.tr("Wavelength of high frequency noise to be filtered\nSet to 4x line spacing")
        )

        self.dlg.lineEdit_3_DC_scale.setToolTip(
            self.tr("Multiplier to be applied to result before subtracting from original grid")
        )

        self.dlg.doubleSpinBox_wtmm_spacing.setToolTip(
            self.tr("Define spacing along profile\n0 = median spacing for points\nfor polylines must be non-zero, and ideally greater than grid cell size")
        )
        self.dlg.pushButton_wtmm.setToolTip(self.tr("Calculate WTMM and display in new windows"))

        self.dlg.mMapLayerComboBox_selectVectors.setToolTip(
            self.tr("Select data points layer or polyline to extract data from grid")
        )

        self.dlg.mFieldComboBox_feature.setToolTip(
            self.tr("Select object to be analysed: LINE_ID of points or FID of polyline")
        )

        self.dlg.mFieldComboBox_data.setToolTip(self.tr("Select data field for points layer"))

        self.dlg.checkBox_PCA.setToolTip(
            self.tr("Principal Component Analysis (PCA)\nOnly works on multiband grids\nReduces dimensionality of data while preserving variance\nUseful for identifying patterns and trends in large datasets")
        )
        self.dlg.mQgsSpinBox_PCA.setToolTip(
            self.tr("Number of components to keep after PCA\nSet to 0 to keep all components")
        )

        self.dlg.checkBox_ICA.setToolTip(
            self.tr("Independent Component Analysis (ICA)\nOnly works on multiband grids\nSeparates a multivariate signal into additive, independent components\nUseful for separating mixed signals and identifying underlying sources")
        )
        self.dlg.mQgsSpinBox_ICA.setToolTip(
            self.tr("Number of components to keep after ICA\nSet to 0 to keep all components")
        )
        self.dlg.checkBox_SunShading.setToolTip(
            self.tr("Sun shading of grid\nUses azimuth and zenith angles to create a shaded relief effect\nUseful for visualizing topography and enhancing features")
        )
        self.dlg.lineEdit_SunSh_Az.setToolTip(
            self.tr("Azimuth angle for sun shading\n0 = North, 90 = East, 180 = South, 270 or -90 = West")
        )
        self.dlg.lineEdit_SunSh_Zn.setToolTip(
            self.tr("Zenith angle for sun shading\n90 = directly overhead, 0 = horizon")
        )
        self.dlg.checkBox_relief.setToolTip(
            self.tr("Uses Grass-like shading algorithm for softer shading")
        )

    def initParams(self):
        """Initialize parameters for the plugin."""
        self.localGridName = ""
        self.diskGridPath = ""
        self.diskPointsPath = ""
        self.buffer = 0
        self.DirClean = False
        self.DC_azimuth = 0
        self.DC_lineSpacing = 400
        self.RTE_P = False
        self.RTE_P_type = "RTP"
        self.RTE_P_inc = 0
        self.RTE_P_dec = 0
        self.RTE_P_height = 0
        self.RTE_P_date = [1, 1, 2000]
        self.RemRegional = False
        self.remReg_wavelength = 5000
        self.Derivative = False
        self.derive_direction = "z"
        self.derive_power = 1.0
        self.TA = False
        self.AS = False
        self.Continuation = False
        self.cont_direction = "up"
        self.cont_height = 500
        self.BandPass = False
        self.band_low = 5
        self.band_high = 50
        self.AGC = False
        self.agc_window = 10
        self.FreqCut = False
        self.FreqCut_type = "Low"
        self.FreqCut_cut = 1000
        self.VI = False
        self.THG = False

        self.Mean = False
        self.mean_size = 3
        self.Median = False
        self.median_size = 3
        self.Gaussian = False
        self.gauss_rad = 1
        self.Direction = "N"
        self.SunShade = False
        self.sun_shade_az = 45
        self.sun_shade_zn = 45

        self.pointType = "point"
        self.input_directory = ""
        self.output_directory = ""
        self.pts_columns = []
        self.pts_data = []

    def parseParams(self):
        """Parse parameters from the dialog."""
        self.DirClean = self.dlg.checkBox_3_DirClean.isChecked()
        self.DC_azimuth = self.dlg.lineEdit_3_azimuth.text()
        self.DC_lineSpacing = self.dlg.lineEdit_3_DC_wavelength.text()
        self.DC_scale = self.to_float(self.dlg.lineEdit_3_DC_scale.text())

        self.RTE_P = self.dlg.checkBox_4_RTE_P.isChecked()
        self.RTE_P_type = self.dlg.comboBox_3_rte_p_list.currentText()
        self.RTE_P_inc = self.to_float(self.dlg.lineEdit_6_inc.text())
        self.RTE_P_dec = self.to_float(self.dlg.lineEdit_5_dec.text())
        self.RTE_P_int = self.to_float(self.dlg.lineEdit_6_int.text())
        date_text = str(self.dlg.dateEdit.date().toPyDate())
        date_split = date_text.split("-")
        self.RTE_P_date = [int(date_split[2]), int(date_split[1]), int(date_split[0])]

        self.RemRegional = self.dlg.checkBox_5_regional.isChecked()
        if self.dlg.radioButton_RR_1st.isChecked():
            self.RemRegional_order = 1
        else:
            self.RemRegional_order = 2

        self.Derivative = self.dlg.checkBox_6_derivative.isChecked()
        self.derive_direction = self.dlg.comboBox_derivDirection.currentText()
        self.derive_power = self.dlg.lineEdit_9_derivePower.text()

        self.TA = self.dlg.checkBox_7_tiltDerivative.isChecked()

        self.AS = self.dlg.checkBox_8_analyticSignal.isChecked()

        self.Continuation = self.dlg.checkBox_9_continuation.isChecked()
        self.cont_direction = self.dlg.comboBox_2_continuationDirection.currentText()
        self.cont_height = self.dlg.lineEdit_10_continuationHeight.text()

        self.BandPass = self.dlg.checkBox_10_bandPass.isChecked()
        self.band_low = self.dlg.lineEdit_12_bandPassLow.text()
        if float(self.band_low) <= 0.0:
            self.band_low = 1e-10
        self.band_high = self.dlg.lineEdit_11_bandPassHigh.text()
        self.band_width = self.dlg.lineEdit_3_BP_width.text()
        self.AGC = self.dlg.checkBox_11_1vd_agc.isChecked()
        self.agc_window = self.dlg.lineEdit_13_agc_window.text()

        self.FreqCut = self.dlg.checkBox_10_freqCut.isChecked()
        self.FreqCut_type = self.dlg.comboBox_2_FreqCutType.currentText()
        self.FreqCut_cut = self.dlg.lineEdit_12_FreqPass.text()
        self.FreqCut_width = self.dlg.lineEdit_3_HLP_width.text()

        self.VI = self.dlg.checkBox_4_PGrav.isChecked()
        self.THG = self.dlg.checkBox_11_tot_hz_grad.isChecked()

        self.Mean = self.dlg.checkBox_Mean.isChecked()
        self.mean_conv_size = self.to_int(self.dlg.lineEdit_Mean_size.text())

        self.Median = self.dlg.checkBox_Median.isChecked()
        self.median_conv_size = self.to_int(self.dlg.lineEdit_Median_size.text())

        self.Gaussian = self.dlg.checkBox_Gaussian.isChecked()
        self.gauss_rad = self.to_float(self.dlg.lineEdit_Gaussian_Sigma.text())

        self.Direction = self.dlg.checkBox_Directional.isChecked()
        self.directional_dir = self.dlg.comboBox_Dir_dir.currentText()

        self.SunShade = self.dlg.checkBox_SunShading.isChecked()
        self.sun_shade_az = self.to_float(self.dlg.lineEdit_SunSh_Az.text())
        self.sun_shade_zn = self.to_float(self.dlg.lineEdit_SunSh_Zn.text())

        self.NaN = self.dlg.checkBox_NaN.isChecked()
        if self.dlg.radioButton_NaN_Above.isChecked():
            self.NaN_Condition = "above"
        elif self.dlg.radioButton_NaN_Below.isChecked():
            self.NaN_Condition = "below"
        else:
            self.NaN_Condition = "between"
        self.NaN_Above = self.to_float(self.dlg.doubleSpinBox_NaN_Above.text())
        self.NaN_Below = self.to_float(self.dlg.doubleSpinBox_NaN_Below.text())

        self.Polygons = self.dlg.checkBox_polygons.isChecked()

        self.SS_Min = self.dlg.checkBox_SS_Min.isChecked()
        self.SS_Max = self.dlg.checkBox_SS_Max.isChecked()
        self.SS_StdDev = self.dlg.checkBox_SS_StdDev.isChecked()
        self.SS_Variance = self.dlg.checkBox_SS_Variance.isChecked()
        self.SS_Skewness = self.dlg.checkBox_SS_Skewness.isChecked()
        self.SS_Kurtosis = self.dlg.checkBox_SS_Kurtosis.isChecked()
        self.SS_Anisotropy = self.dlg.checkBox_SS_Anisotropy.isChecked()
        win_type_text = self.dlg.comboBox_SS_Anisotropy_WinType.currentText()
        self.SS_anisotropy_window_type = 'gaussian' if win_type_text == self.tr("Gaussian") else 'box'
        self.SS_ChainLength = self.dlg.checkBox_SS_ChainLength.isChecked()
        self.SS_Streamline  = self.dlg.checkBox_SS_Streamline.isChecked()
        self.SS_aniso_threshold = float(self.dlg.doubleSpinBox_SS_AnisoThresh.value())
        self.SS_angle_tolerance = float(self.dlg.doubleSpinBox_SS_AngleTol.value())
        self.SS_max_steps       = int(self.dlg.spinBox_SS_MaxSteps.value())
        self.SS_search_radius   = int(self.dlg.spinBox_SS_SearchRadius.value())
        self.SS_window_size = self.to_int(self.dlg.lineEdit_SS_Window.text())

        self.MRVBF = self.dlg.checkBox_MRVBF.isChecked()
        self.MRVBF_t_slope = self.to_float(self.dlg.lineEdit_MRVBF_tslope.text())
        self.MRVBF_t_pctl_v = self.to_float(self.dlg.lineEdit_MRVBF_tpctlv.text())
        self.MRVBF_t_pctl_r = self.to_float(self.dlg.lineEdit_MRVBF_tpctlr.text())
        self.MRVBF_p_slope = self.to_float(self.dlg.lineEdit_MRVBF_pslope.text())
        self.MRVBF_p_pctl = self.to_float(self.dlg.lineEdit_MRVBF_ppctl.text())
        self.MRVBF_max_res = self.to_float(self.dlg.lineEdit_MRVBF_maxres.text())

        self.PCA = self.dlg.checkBox_PCA.isChecked()
        self.ICA = self.dlg.checkBox_ICA.isChecked()
        self.ED = self.dlg.checkBox_ED.isChecked()

    def loadGrid(self):
        """
        Loads a raster grid from the specified file path, adds it to the QGIS project if not already loaded,
        and processes the raster data into a NumPy array.
        This method performs the following steps:
        1. Loads the raster file as a QgsRasterLayer.
        2. Adds the layer to the QGIS project if it is not already loaded.
        3. Retrieves raster metadata such as pixel dimensions and raster extent.
        4. Reads the raster data into a NumPy array.
        5. Handles NoData values by replacing them with NaN in the NumPy array.
        Attributes:
            self.layer (QgsRasterLayer): The loaded raster layer.
            self.dx (float): The pixel width in map units.
            self.dy (float): The pixel height in map units.
            self.raster_array (numpy.ndarray): The raster data as a 2D NumPy array.
        Raises:
            ValueError: If the raster data cannot be loaded or processed.
        Notes:
            - The method assumes that `self.diskGridPath` contains the valid file path to the raster file.
            - The raster data is read from the first band (band 1).
            - NoData values are replaced with NaN in the resulting NumPy array.
        """

        fileInfo = QFileInfo(self.diskGridPath)
        baseName = fileInfo.baseName()

        self.layer = QgsRasterLayer(self.diskGridPath, baseName)
        if not self.is_layer_loaded(baseName):
            QgsProject.instance().addMapLayer(self.layer)

        self.dx = self.layer.rasterUnitsPerPixelX()
        self.dy = self.layer.rasterUnitsPerPixelY()
        """
        # Access the raster data provider
        provider = self.layer.dataProvider()

        # Get raster dimensions
        cols = provider.xSize()  # Number of columns
        rows = provider.ySize()  # Number of rows

        # Read raster data as a block
        band = 1  # Specify the band number (1-based index)
        raster_block = provider.block(band, provider.extent(), cols, rows)

        # Copy the block data into a NumPy array
        extent = self.layer.extent()
        rows, cols = self.layer.height(), self.layer.width()
        raster_block = provider.block(1, extent, cols, rows)  # !!!!!
        self.raster_array = np.zeros((rows, cols))
        for i in range(rows):
            for j in range(cols):
                self.raster_array[i, j] = raster_block.value(i, j)
        """
        # Read raster data via GDAL (consistent row ordering across platforms)
        ds = gdal.Open(self.diskGridPath)
        band = ds.GetRasterBand(1)
        no_data_value = band.GetNoDataValue()  # Band 1
        nodata = ds.GetRasterBand(1).GetNoDataValue()
        band_data = ds.GetRasterBand(1).ReadAsArray()  # Always row0=north
        if nodata is not None:
            band_data = band_data.astype(float)
            band_data[band_data == nodata] = np.nan
        self.raster_array = band_data
        gt = ds.GetGeoTransform()
        if gt[5] > 0:  # positive pixel height means south-up, flip to north-up
            self.raster_array = np.flipud(self.raster_array)
        ds = None

        # Handle NoData values if needed

        if no_data_value is not None:
            self.raster_array[self.raster_array == no_data_value] = np.nan

    def insert_text_before_extension(self, file_path, insert_text):
        """
        Insert text at the end of the filename, before the file extension.

        Parameters:
            file_path (str): Full path of the file.
            insert_text (str): Text to insert before the file extension.

        Returns:
            str: The modified file path.
        """
        # Separate the file path into directory, base name, and extension
        dir_name, base_name = os.path.split(file_path)
        file_name, file_ext = os.path.splitext(base_name)

        # Construct the new file name
        new_file_name = f"{file_name}{insert_text}{file_ext}"

        # Combine directory and new file name
        return os.path.join(dir_name, new_file_name)

    def procIDWGridding(self):
        gridder = QGISGridData(self.iface)

        layer_name = self.dlg.mMapLayerComboBox_selectGrid_3.currentText()
        input = self.get_layer_path_by_name(layer_name)
        zcolumn = self.dlg.comboBox_select_grid_data_field.currentText()
        cell_size = self.dlg.doubleSpinBox_cellsize.text()

        mask = None
        alg_id = "grass7:v.surf.idw"
        try:
            # Check if the algorithm exists
            if QgsApplication.processingRegistry().algorithmById(alg_id):
                # Launch the dialog
                gridder.launch_idw_dialog(input, zcolumn, cell_size, mask)
            else:
                self.iface.messageBar().pushMessage(
                    "Error",
                    "GRASS v.surf.rst algorithm not found.",
                    level=Qgis.MessageLevel.Critical,
                )
        except Exception as e:
            self.iface.messageBar().pushMessage(
                "Error", str(e), level=Qgis.MessageLevel.Critical
            )

    @staticmethod
    def _mba_grid_geometry(data_xmin, data_xmax, data_ymin, data_ymax, cell_size):
        """SAGA Multilevel B-Spline grid target: the lower-left node is snapped
        DOWN to a multiple of the cell size and the node count is
        2 + floor(dataRange / cellSize) -- one node for the snapped minimum,
        one per whole cell of data range, plus one to span the partial last
        cell. Returns (xmin, ymin, nx, ny) where xmin/ymin are the coordinates
        of cell (0, 0)'s centre. Used both for the actual gridding and for the
        '# Cells' preview labels so they agree with the B-spline output."""
        xmin = float(np.floor(data_xmin / cell_size) * cell_size)
        ymin = float(np.floor(data_ymin / cell_size) * cell_size)
        nx = 2 + int((data_xmax - xmin) / cell_size)
        ny = 2 + int((data_ymax - ymin) / cell_size)
        return xmin, ymin, nx, ny

    def procmultibsplineGridding(self):
        """Multilevel B-Spline (MBA) gridding of the selected point layer /
        data field into a temporary GeoTIFF, using the local NumPy translation
        of SAGA's grid_spline tool (calcs/saga_mba_gridding.py).

        Grid geometry (extent, cell size) and the value field are taken from
        the Grid + Wavelets tab widgets. The result is loaded into the project.
        """
        layer_name = self.dlg.mMapLayerComboBox_selectGrid_3.currentText()
        if not layer_name:
            self.iface.messageBar().pushMessage(
                "Error", self.tr("No points layer selected"),
                level=Qgis.MessageLevel.Critical,
            )
            return

        layers = QgsProject.instance().mapLayersByName(layer_name)
        if not layers:
            return
        layer = layers[0]

        zcolumn = self.dlg.comboBox_select_grid_data_field.currentText()
        if not zcolumn or zcolumn not in [f.name() for f in layer.fields()]:
            self.iface.messageBar().pushMessage(
                "Error", self.tr("Select a valid data field to grid"),
                level=Qgis.MessageLevel.Critical,
            )
            return

        cell_size = float(self.dlg.doubleSpinBox_cellsize.value())
        epsilon = self.to_float(self.dlg.lineEdit_bspline_epsilon.text())
        level_max = int(self.dlg.spinBox_bspline_levels.value())

        # ---- collect scattered points ----------------------------------
        xs, ys, zs = [], [], []
        for feat in layer.getFeatures():
            geom = feat.geometry()
            if geom is None or geom.isEmpty():
                continue
            val = feat[zcolumn]
            if val is None:
                continue
            try:
                zval = float(val)
            except (TypeError, ValueError):
                continue
            if geom.isMultipart():
                pts = geom.asMultiPoint()
            else:
                pts = [geom.asPoint()]
            for pt in pts:
                xs.append(pt.x())
                ys.append(pt.y())
                zs.append(zval)

        if len(xs) < 3:
            self.iface.messageBar().pushMessage(
                "Error",
                self.tr("Need at least 3 valid points with a numeric value to grid"),
                level=Qgis.MessageLevel.Critical,
            )
            return

        xs = np.asarray(xs, dtype=float)
        ys = np.asarray(ys, dtype=float)
        zs = np.asarray(zs, dtype=float)

        # ---- optional nodata-sentinel filter -------------------------
        if self.dlg.checkBox_bspline_ignore.isChecked():
            thresh = self.to_float(self.dlg.lineEdit_bspline_ignore.text())
            keep = zs >= thresh
            dropped = int((~keep).sum())
            xs, ys, zs = xs[keep], ys[keep], zs[keep]
            if dropped:
                self.iface.messageBar().pushMessage(
                    "B-Spline gridding: ignored %d point(s) with %s < %g"
                    % (dropped, zcolumn, thresh),
                    level=Qgis.Info,
                    duration=8,
                )
            if xs.size < 3:
                self.iface.messageBar().pushMessage(
                    "Error",
                    self.tr("Too few points remain after the ignore-value filter"),
                    level=Qgis.MessageLevel.Critical,
                )
                return

        # ---- detrend by the data minimum -----------------------------
        # Geophysical fields often sit on a large offset with a small dynamic
        # range (e.g. ~57700 nT +/- 100). Grid the residual and restore the
        # offset afterwards so the B-spline solve keeps its precision.
        # (mba_gridding also mean-centres internally; this is belt-and-braces
        # and makes the behaviour explicit.)
        # z_offset = float(np.min(zs))
        # zs = zs - z_offset

        # ---- target grid geometry -----------------------------------
        xmin, ymin, nx, ny = self._mba_grid_geometry(
            float(xs.min()), float(xs.max()),
            float(ys.min()), float(ys.max()), cell_size,
        )
        if nx < 2 or ny < 2:
            self.iface.messageBar().pushMessage(
                "Error", self.tr("Cell size too large for the point extent"),
                level=Qgis.MessageLevel.Critical,
            )
            return

        try:
            grid = mba_gridding(
                xs, ys, zs,
                cellsize=cell_size,
                xmin=xmin, ymin=ymin,
                nx=nx, ny=ny,
                epsilon=epsilon,
                level_max=level_max,
                refinement=False,
            )
        except Exception as e:
            self.iface.messageBar().pushMessage(
                "Error", str(e), level=Qgis.MessageLevel.Critical
            )
            return

        # grid = grid + z_offset  # restore the removed offset

        # mba_gridding returns row 0 = ymin; GeoTIFF is written north-up.
        grid = np.flipud(grid)

        out_path = os.path.join(
            tempfile.gettempdir(),
            f"{layer_name}_{zcolumn}_bspline.tif",
        )
        if os.path.exists(out_path):
            try:
                os.remove(out_path)
            except OSError:
                out_path = tempfile.NamedTemporaryFile(
                    suffix=".tif", delete=False
                ).name

        rows, cols = grid.shape
        driver = gdal.GetDriverByName("GTiff")
        ds = driver.Create(out_path, cols, rows, 1, gdal.GDT_Float32)
        # mba_gridding places cell (0,0) *centre* at (xmin, ymin); GDAL's
        # geotransform origin is the top-left pixel *corner*.
        gt_originx = xmin - cell_size / 2.0
        gt_originy = ymin + (ny - 0.5) * cell_size
        ds.SetGeoTransform([gt_originx, cell_size, 0, gt_originy, 0, -cell_size])
        srs = osr.SpatialReference()
        authid = layer.crs().authid()
        if authid and ":" in authid:
            srs.ImportFromEPSG(int(authid.split(":")[1]))
            ds.SetProjection(srs.ExportToWkt())
        band = ds.GetRasterBand(1)
        band.SetNoDataValue(np.nan)
        band.WriteArray(grid.astype(np.float32))
        band.FlushCache()
        band = None
        ds = None

        layer_out_name = f"{layer_name}_{zcolumn}_bspline"
        if self.is_layer_loaded(layer_out_name):
            project = QgsProject.instance()
            for lyr in project.mapLayersByName(layer_out_name):
                project.removeMapLayer(lyr.id())

        raster_layer = QgsRasterLayer(out_path, layer_out_name)
        if raster_layer.isValid():
            QgsProject.instance().addMapLayer(raster_layer)
        else:
            self.iface.messageBar().pushMessage(
                "Error", self.tr("Failed to load gridded result"),
                level=Qgis.MessageLevel.Critical,
            )

    def get_layer_path_by_name(self, layer_name):
        """
        Get the file path of a layer given its name.

        :param layer_name: The name of the layer in the QGIS project.
        :return: File path of the layer or None if not found or not a file-based layer.
        """
        # Iterate through all layers in the project
        for layer in QgsProject.instance().mapLayers().values():
            if layer.name() == layer_name:
                # Check if the layer has a data provider and a file source
                if hasattr(layer, "dataProvider") and hasattr(
                    layer.dataProvider(), "dataSourceUri"
                ):
                    # Return the file path (for file-based layers like shapefiles or rasters)
                    return layer.dataProvider().dataSourceUri().split("|")[0]
        return None

    def procDirClean(self):
        cutoff_wavelength = 4 * float(self.DC_lineSpacing)
        # if self.unit_check(cutoff_wavelength) or True:
        self.new_grid = self.processor.directional_butterworth_band_pass(
            self.raster_array,
            1e-8,
            float(self.DC_lineSpacing),
            direction_angle=float(self.DC_azimuth),
            direction_width=20,
            order=4,
            buffer_size=10,
            buffer_method="mirror",
        )
        print("xxxxx")
        nan_mask = np.isnan(self.new_grid)
        self.new_grid[nan_mask] = 1.0
        self.new_grid = self.raster_array - (self.new_grid * self.DC_scale)
        self.new_grid[nan_mask] = np.nan
        self.suffix = "_DirC"

    def procRTP_E(self, sgtool_instance):
        if self.RTE_P_inc == "0" and self.RTE_P_dec == "0":
            self.iface.messageBar().pushMessage(
                "You need to define Inc and Dec first!", level=Qgis.Warning, duration=15
            )
        else:
            if self.RTE_P_type == "RTP":
                self.new_grid = self.processor.reduction_to_pole(
                    self.raster_array,
                    inclination=float(self.RTE_P_inc),
                    declination=float(self.RTE_P_dec),
                    buffer_size=self.buffer,
                )
                self.suffix = "_RTP"
            else:
                self.new_grid = self.processor.reduction_to_equator(
                    self.raster_array,
                    inclination=float(self.RTE_P_inc),
                    declination=float(self.RTE_P_dec),
                    buffer_size=self.buffer,
                )
                self.new_grid = self.new_grid  # * float(self.RTE_P_int)
                self.suffix = "_RTE"

    def procVRTP(self):
        """Differential (variable) RTP using Cooper/Cowan Taylor-series method.

        Uses the Cooper/Cowan rtp() for all 13 Taylor-series calls,
        giving results identical to the reference standalone implementation.
        """
        inc_corners, dec_corners, inc_center, dec_center = self.get_igrf_corners()
        if inc_corners is None:
            self.iface.messageBar().pushMessage(
                "Var. RTP: could not compute IGRF corners — check layer CRS and date.",
                level=Qgis.Warning, duration=15,
            )
            return

        data = self.raster_array
        nr, nc = data.shape

        # SGTool stores rasters north-up (row 0 = north). Cooper/Cowan expects
        # MATLAB axis-xy convention (row 1 = south, row nr = north). Flip before
        # passing in so incv is interpolated in the correct geographic direction,
        # then flip the result back to north-up for QGIS display.
        # NaN filling and restoration are handled internally by rtpvariable.
        _, varrtp, incv_deg, decv_deg = cooper_rtpvariable(
            np.flipud(data), nr, nc, inc_corners, dec_corners,
            inc_center=inc_center, dec_center=dec_center,
        )
        self.new_grid = np.flipud(varrtp)
        self.suffix = "_DRTP"

        # Save inc/dec field grids to disk alongside the output (not added to QGIS)
        # Flip south-up grids back to north-up to match GeoTIFF convention
        for grid, tag in ((np.flipud(incv_deg), "_DRTP_inc"), (np.flipud(decv_deg), "_DRTP_dec")):
            out_path = self.insert_text_before_extension(self.diskGridPath, tag)
            if ".tif" not in out_path.lower():
                out_path = os.path.splitext(out_path)[0] + ".tif"
            self.numpy_array_to_raster(grid, out_path, reference_layer=self.layer)

    def procRemRegional(self):

        data, nodata_value = self.processor.fix_extreme_values(self.raster_array)
        data = data.astype(np.float32)
        mask = (data == nodata_value) | np.isnan(data)

        if self.RemRegional_order == 1:
            self.new_grid = self.processor.remove_gradient(data, mask)
        else:
            self.new_grid = self.processor.remove_2o_gradient(data, mask)

        self.suffix = "_RR" + "_" + str(self.RemRegional_order) + "o"

    def procDerivative(self):
        self.new_grid = self.processor.compute_derivative(
            self.raster_array,
            direction=self.derive_direction,
            order=float(self.derive_power),
            buffer_size=self.buffer,
        )
        self.suffix = "_d" + str(self.derive_power) + self.derive_direction

    def procTiltAngle(self):
        self.new_grid = self.processor.tilt_angle(
            self.raster_array, buffer_size=self.buffer
        )
        self.suffix = "_TA"

    def procAnalyticSignal(self):
        self.new_grid = self.processor.analytic_signal(
            self.raster_array, buffer_size=self.buffer
        )
        self.suffix = "_AS"

    def procContinuation(self):
        selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]

        crs = selected_layer.crs()
        if crs.isGeographic():
            long, lat = self.get_grid_centroid(selected_layer)
            dx, dy = self.SG_Util.arc_degree_to_meters(lat)
            ave_dxdy = np.sqrt(dx**2.0 + dy**2.0) / 2

            height = float(self.cont_height) / ave_dxdy
            self.iface.messageBar().pushMessage(
                "Height roughly converted to metres, since this is a geographic projection",
                level=Qgis.Success,
                duration=15,
            )
        else:
            height = float(self.cont_height)
        if self.cont_direction == "up":
            self.new_grid = self.processor.upward_continuation(
                self.raster_array, height=height, buffer_size=self.buffer
            )
            self.suffix = "_UC" + "_" + str(self.cont_height)
        else:
            self.new_grid = self.processor.downward_continuation(
                self.raster_array, height=height, buffer_size=self.buffer
            )
            self.suffix = "_DC" + "_" + str(self.cont_height)

    def procBandPass(self):
        low_cut = float(self.band_low)
        high_cut = float(self.band_high)
        if self.unit_check(low_cut) and self.unit_check(high_cut):
            """self.new_grid = self.processor.butterworth_band_pass(
                self.raster_array, low_cut, high_cut, order=4, buffer_size=self.buffer
            )"""
            self.new_grid = self.processor.band_pass_filter(
                self.raster_array,
                low_cut=low_cut,
                high_cut=high_cut,
                high_transition_width=float(self.band_width),
                low_transition_width=float(self.band_width),
                buffer_size=self.buffer,
            )
            self.suffix = "_BP" + "_" + str(self.band_low) + "_" + str(self.band_high)

    def procAGC(self):
        self.new_grid = self.processor.automatic_gain_control(
            self.raster_array, window_size=int(self.agc_window)
        )
        self.suffix = "_AGC"

    def procTHG(self):
        self.new_grid = self.processor.total_hz_grad(
            self.raster_array, buffer_size=self.buffer
        )
        self.suffix = "_THG"

    def procvInt(self):
        selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
        crs = selected_layer.crs()
        if crs.isGeographic():
            self.iface.messageBar().pushMessage(
                "Vertical integration requires a metre-based projection system",
                level=Qgis.Warning,
                duration=15,
            )
        else:
            self.new_grid = self.processor.vertical_integration(
                self.raster_array,
                max_wavenumber=None,
                min_wavenumber=1e-4,
                buffer_size=self.buffer,
                buffer_method="mirror",
            )
            self.suffix = "_VI"

    def procFreqCut(self):
        cutoff_wavelength = float(self.FreqCut_cut)
        if self.unit_check(cutoff_wavelength):
            if self.FreqCut_type == "Low":
                self.new_grid = self.processor.low_pass_filter(
                    self.raster_array,
                    cutoff_wavelength=cutoff_wavelength,
                    transition_width=float(self.FreqCut_width),
                    buffer_size=self.buffer,
                )
                self.suffix = "_LP" + "_" + str(self.FreqCut_cut)
            else:
                self.new_grid = self.processor.high_pass_filter(
                    self.raster_array,
                    cutoff_wavelength=cutoff_wavelength,
                    transition_width=float(self.FreqCut_width),
                    buffer_size=self.buffer,
                )
                self.suffix = "_HP" + "_" + str(self.FreqCut_cut)

    def procMean(self):
        self.new_grid = self.convolution.mean_filter(
            # self.raster_array,
            self.mean_conv_size
        )
        self.suffix = "_Mn"

    def procMedian(self):
        self.new_grid = self.convolution.median_filter(
            # self.raster_array,
            self.median_conv_size
        )
        self.suffix = "_Md"

    def procGaussian(self):
        self.new_grid = self.convolution.gaussian_filter(
            # self.raster_array,
            self.gauss_rad
        )
        self.suffix = "_Gs"

    def procDirectional(self):
        self.new_grid = self.convolution.directional_filter(
            # self.raster_array,
            self.directional_dir,
            n=3,
        )
        self.suffix = "_Dr"

    def procSunShade(self):

        if self.dlg.checkBox_relief.isChecked():
            selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[
                0
            ]
            crs = selected_layer.crs()
            if crs.isGeographic():
                long, lat = self.get_grid_centroid(selected_layer)
                dx, dy = self.SG_Util.arc_degree_to_meters(lat)
                ave_dxdy = np.sqrt(dx**2.0 + dy**2.0) / 2
                hzscale = 1 * ave_dxdy
            else:
                hzscale = 1.0

            self.new_grid = self.convolution.sun_shading_filter_grass(
                self.raster_array,
                altitude=self.sun_shade_zn,
                azimuth=self.sun_shade_az,
                resolution_ns=self.dy * hzscale,
                resolution_ew=self.dx * hzscale,
                scale=1.0,
                zscale=1.0,
            )
        else:
            self.new_grid = self.convolution.sun_shading_filter(
                self.raster_array,
                sun_alt=self.sun_shade_zn,
                sun_az=180 - self.sun_shade_az,
            )
        self.suffix = "_Sh"

    def procEulerDeconvolution(self, data):
        selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
        self.diskGridPath = selected_layer.dataProvider().dataSourceUri()
        crs = selected_layer.crs()

        if crs.isGeographic():
            self.iface.messageBar().pushMessage(
                "This is a geographic projection, you need to convert it to a projected CRS",
                level=Qgis.Warning,
                duration=15,
            )
            return False
        else:
            print("Processing NaNs")
            data, mask = self.processor.fill_nan(data)

            shape = (data.shape[0], data.shape[1])
            provider = selected_layer.dataProvider()

            # Get raster dimensions
            area = provider.extent()
            # use south, north, west, east order
            area = [area.yMinimum(), area.yMaximum(), area.xMinimum(), area.xMaximum()]
            # print("area", area)
            # Assuming your data array is called 'data'
            rows, cols = data.shape
            # print("rows,cols", rows, cols)
            # print("data.shape", data.shape)
            # Create 1D coordinate arrays
            x_coords_1d = np.linspace(area[2], area[3], cols)
            y_coords_1d = np.linspace(area[0], area[1], rows)

            # Create 2D coordinate grids
            yi, xi = np.meshgrid(x_coords_1d, y_coords_1d)
            zi = np.ones(shape)
            # zi = zi * -100.0
            xi = xi.flatten()
            yi = yi.flatten()
            zi = zi.flatten()
            data = np.flipud(data)
            data = data.flatten()
            # print("data,xi,yi,zi", data.shape, xi.shape, yi.shape, zi.shape)
            # print("data[:200]", data[:200])
            # print("xi[:200]", xi[:200])
            # print("yi[:200]", yi[:200])
            # print("zi[:200]", zi[:200])
            # moving data window size
            winsize = self.to_int(self.dlg.lineEdit_ED_Window.text())
            # percentage of the solutions that will be keep
            filt = self.to_float(self.dlg.doubleSpinBox_ED_Threshold.text())

            # print("winsize,filt", winsize, filt)
            # empty array for multiple SIs
            est_classic = []
            # Define below the SIs to be tested
            SI_vet = [0.001, 1, 2, 3]
            # prism (SI = 0), a line of poles (SI = 1), a single pole (SI = 2), and a di-pole (SI = 3)
            """print(
                "data, xi, yi, zi, shape, area, SI, winsize, filt",
                data.shape,
                xi.shape,
                yi.shape,
                zi.shape,
                shape,
                area,
                SI_vet,
                winsize,
                filt,
            )"""
            """
            Euler deconvolution for multiple SIs
            """
            for SI in SI_vet:
                print(f"Processing Euler Deconvolution for SI = {SI}")
                classic_result = euler_deconv(
                    data, xi, yi, zi, shape, area, SI, winsize, filt
                )
                classic_result[1, :] = (
                    area[0] - classic_result[1, :] + area[1]
                )  # Flip x-coordinates
                est_classic.append(classic_result)

            area_classic = area

            head_tail = os.path.split(self.diskGridPath)
            # convert the list to an array
            output = np.asarray(est_classic)
            # save the estimates in distinct files according to the SI
            for i in range(len(SI_vet)):
                np.savetxt(
                    head_tail[0]
                    + "/"
                    + self.localGridName
                    + "_estimates_SI_"
                    + str(i)
                    + ".txt",
                    output[i],
                    delimiter=",",
                    header="y_source, x_source, z_source, base_level",
                    comments="",  # This removes the # prefix
                )
            # optional windowed stats
            if self.dlg.checkBox_ED_Stats.isChecked():
                # classic(est_classic, area_classic, SI_vet, self.localGridName, head_tail[0])
                window_results = window_stats(
                    est_classic,
                    area,
                    SI_vet,
                    self.localGridName,
                    head_tail[0],
                    shape,
                    winsize,
                    detailed_stats=True,
                )

            self.iface.messageBar().pushMessage(
                "Euler Solutions saved to same directory as input grid",
                level=Qgis.Success,
                duration=15,
            )

    def procPCA(self):
        try:
            import sklearn
        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for PCAs, but not installed: scikit-learn\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install scikit-learn",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False
        self.suffix = "_PCA"
        selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
        self.diskGridPath = selected_layer.dataProvider().dataSourceUri()
        self.diskNewGridPath = self.insert_text_before_extension(
            self.diskGridPath, self.suffix
        )
        n_components = int(self.dlg.mQgsSpinBox_PCA.value())
        self.delete_layer_and_file(self.localGridName + self.suffix)

        raster_layer = QgsRasterLayer(self.diskGridPath, "Input Raster")
        if not raster_layer.isValid():
            raise ValueError(f"Failed to load raster layer: {self.diskGridPath}")

        components, variance_ratio = self.PCAICA.pca_with_nans(
            self.diskGridPath, self.diskNewGridPath, n_components
        )
        if components is not None:
            PCA_raster_layer = QgsRasterLayer(
                self.diskNewGridPath, self.localGridName + self.suffix
            )
            QgsProject.instance().addMapLayer(PCA_raster_layer)

    def delete_layer_and_file(self, layer_name, file_path=None):
        """
        Check if a layer with the given name exists in QGIS and delete it.
        Also check if the file exists on disk and delete it.

        Parameters:
        layer_name (str): The name of the layer in QGIS
        file_path (str, optional): The file path to check and delete. If None,
                                will try to get the path from the layer

        Returns:
        tuple: (layer_deleted, file_deleted) indicating if layer and file were deleted
        """
        project = QgsProject.instance()
        layer_deleted = False
        file_deleted = False

        # Check if layer exists in the project
        layer = project.mapLayersByName(layer_name)
        if layer:
            # Layer exists in the project, get the first match
            layer = layer[0]

            # If file_path wasn't provided, try to get it from the layer
            if file_path is None:
                if hasattr(layer, "source"):
                    file_path = layer.source()

            # Remove layer from the project
            layer_id = layer.id()
            success = project.removeMapLayer(layer_id)
            if success is None:  # removeMapLayer returns None on success
                print(f"Layer '{layer_name}' has been removed from the project")
                layer_deleted = True
            else:
                print(f"Failed to remove layer '{layer_name}' from the project")
        else:
            print(f"Layer '{layer_name}' not found in the project")

        # Check if file exists on disk and delete it
        if file_path:
            if os.path.exists(file_path):
                try:
                    # Check if it's a directory
                    if os.path.isdir(file_path):
                        import shutil

                        shutil.rmtree(file_path)
                        # print(f"Directory '{file_path}' has been deleted")
                        file_deleted = True
                    else:
                        # It's a file
                        os.remove(file_path)
                        # print(f"File '{file_path}' has been deleted")
                        file_deleted = True
                except Exception as e:
                    print(f"Failed to delete '{file_path}': {str(e)}")
            else:
                print(f"File '{file_path}' not found on disk")

        return layer_deleted, file_deleted

    def procICA(self):
        try:
            import sklearn
        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for ICAs, but not installed: scikit-learn\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install scikit-learn",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False
        self.suffix = "_ICA"
        selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]

        self.diskGridPath = selected_layer.dataProvider().dataSourceUri()
        self.diskNewGridPath = self.insert_text_before_extension(
            self.diskGridPath, self.suffix
        )
        n_components = int(self.dlg.mQgsSpinBox_ICA.value())
        self.delete_layer_and_file(self.localGridName + self.suffix)

        raster_layer = QgsRasterLayer(self.diskGridPath, "Input Raster")
        if not raster_layer.isValid():
            raise ValueError(f"Failed to load raster layer: {self.diskGridPath}")

        mixing_matrix, unmixing_matrix = self.PCAICA.ica_with_nans(
            self.diskGridPath, self.diskNewGridPath, n_components
        )
        if mixing_matrix is not None:
            ICA_raster_layer = QgsRasterLayer(
                self.diskNewGridPath, self.localGridName + self.suffix
            )
            QgsProject.instance().addMapLayer(ICA_raster_layer)

    def procPolygons(self):
        if self.localGridName and self.localGridName != "":
            self.parseParams()

            self.layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
            input_raster_path = self.layer.dataProvider().dataSourceUri()
            output_path = self.insert_text_before_extension(
                input_raster_path, "_boundary"
            )
            output_path_shp = self.SG_Util.create_data_boundary_lines(
                input_raster_path, output_path
            )
            layer = QgsVectorLayer(output_path_shp, self.localGridName + "_boundary")
            if not layer.isValid():
                raise ValueError(f"Failed to load layer: {output_path_shp}")
            else:
                # Add the layer to the current QGIS project
                QgsProject.instance().addMapLayer(layer)

    def procNaN(self):
        self.new_grid = self.SG_Util.Threshold2Nan(
            self.raster_array,
            condition=self.NaN_Condition,
            above_threshold_value=self.NaN_Above,
            below_threshold_value=self.NaN_Below,
        )
        self.suffix = "_Clean"

    def procNormalise(self):
        processor = GeophysicalProcessor(None, None, None)
        inpath = self.input_directory
        outpath = self.output_directory
        order = self.dlg.radioButton_normalise_1st.isChecked()
        if (
            os.path.exists(inpath)
            and os.path.exists(outpath)
            and inpath != ""
            and outpath != ""
        ):
            processor.normalise_geotiffs(inpath, outpath, order)

    def procSS_Min(self):
        self.new_grid = self.SpatialStats.calculate_windowed_stats(
            window_size=self.SS_window_size, stat_type="min"
        )
        self.suffix = "_SS_Min"

    def procSS_Max(self):
        self.new_grid = self.SpatialStats.calculate_windowed_stats(
            window_size=self.SS_window_size, stat_type="max"
        )
        self.suffix = "_SS_Max"

    def procSS_StdDev(self):
        self.new_grid = self.SpatialStats.calculate_windowed_stats(
            window_size=self.SS_window_size, stat_type="std"
        )
        self.suffix = "_SS_StdDev"

    def procSS_Variance(self):
        self.new_grid = self.SpatialStats.calculate_windowed_stats(
            window_size=self.SS_window_size, stat_type="variance"
        )
        self.suffix = "_SS_Var"

    def procSS_Skewness(self):
        self.new_grid = self.SpatialStats.calculate_windowed_stats(
            window_size=self.SS_window_size, stat_type="skewness"
        )
        self.suffix = "_SS_Skew"

    def procSS_Kurtosis(self):
        self.new_grid = self.SpatialStats.calculate_windowed_stats(
            window_size=self.SS_window_size, stat_type="kurtosis"
        )
        self.suffix = "_SS_Kurt"

    def procSS_Anisotropy(self):
        anisotropy_map, orientation_deg = self.SpatialStats.calculate_local_anisotropy(
            self.SpatialStats.grid,
            window_size=self.SS_window_size,
            window_type=self.SS_anisotropy_window_type,
        )
        self._aniso_orientation = orientation_deg  # stashed for main loop
        self.new_grid = anisotropy_map
        self.suffix = "_SS_AnisoMag"

    def _get_anisotropy_maps(self):
        """Compute anisotropy magnitude and orientation using current settings."""
        return self.SpatialStats.calculate_local_anisotropy(
            self.SpatialStats.grid,
            window_size=self.SS_window_size,
            window_type=self.SS_anisotropy_window_type,
        )

    def procSS_ChainLength(self):
        aniso_map, orient_map = self._get_anisotropy_maps()
        self.new_grid = self.SpatialStats.calculate_anisotropy_chain_length(
            aniso_map, orient_map,
            aniso_threshold=self.SS_aniso_threshold,
            angle_tolerance=self.SS_angle_tolerance,
            search_radius=self.SS_search_radius,
        )
        self.suffix = "_SS_ChainLen"
        self.addNewGrid(stdClip=False)

    def procSS_Streamline(self):
        aniso_map, orient_map = self._get_anisotropy_maps()
        self.new_grid = self.SpatialStats.calculate_anisotropy_streamline_length(
            aniso_map, orient_map,
            aniso_threshold=self.SS_aniso_threshold,
            angle_tolerance=self.SS_angle_tolerance,
            max_steps=self.SS_max_steps,
        )
        self.suffix = "_SS_StreamLen"
        self.addNewGrid(stdClip=False)

    def procMRVBF(self):
        """Multiresolution Valley Bottom Flatness / Ridge Top Flatness.

        Runs the Gallant & Dowling (2003) MRVBF/MRRTF translation on the
        selected grid, computes the slope of the same grid with the built-in
        QGIS slope algorithm, and writes:

            <name>_MRVBF.tif        valley bottom flatness index
            <name>_MRRTF.tif        ridge top flatness index
            <name>_MRVBF_slope.tif  slope (degrees) of the input grid
            <name>_MRVBF_RGB.tif    3-band composite (R=MRRTF, G=MRVBF, B=slope)

        All four are loaded into the project. Geographic (lat/long) rasters
        are rejected because the algorithm needs cell sizes in metres.
        """
        crs = self.layer.crs()
        if crs.isGeographic():
            self.iface.messageBar().pushMessage(
                self.tr(
                    "MRVBF needs a projected raster with cell sizes in metres. "
                    "Reproject the grid to a projected CRS and try again."
                ),
                level=Qgis.Warning,
                duration=15,
            )
            return

        # SGTool keeps the array north-up (row 0 = north); mrvbf() expects
        # row 0 = ymin, so flip going in and flip every result back.
        dem_su = np.flipud(self.raster_array.astype(float))
        nodata_mask = ~np.isfinite(dem_su)

        extent = self.layer.dataProvider().extent()
        cellsize = float(self.dx)

        mrvbf_su, mrrtf_su = calc_mrvbf(
            dem_su,
            cellsize,
            xmin=extent.xMinimum(),
            ymin=extent.yMinimum(),
            nodata_mask=nodata_mask,
            t_slope=self.MRVBF_t_slope,
            t_pctl_v=self.MRVBF_t_pctl_v,
            t_pctl_r=self.MRVBF_t_pctl_r,
            p_slope=self.MRVBF_p_slope,
            p_pctl=self.MRVBF_p_pctl,
            max_res=self.MRVBF_max_res,
        )
        mrvbf_grid = np.flipud(mrvbf_su)
        mrrtf_grid = np.flipud(mrrtf_su)

        # --- MRVBF / MRRTF grids (written + loaded via addNewGrid) ---
        self.new_grid = mrvbf_grid
        self.suffix = "_MRVBF"
        self.addNewGrid(stdClip=False)

        self.new_grid = mrrtf_grid
        self.suffix = "_MRRTF"
        self.addNewGrid(stdClip=False)

        # --- slope of the input grid via the built-in QGIS slope algorithm ---
        slope_path = self.insert_text_before_extension(
            self.diskGridPath, "_MRVBF_slope"
        )
        base, _ext = os.path.splitext(slope_path)
        slope_path = base + ".tif"
        if os.path.exists(slope_path):
            try:
                os.remove(slope_path)
            except OSError:
                pass
        try:
            processing.run(
                "native:slope",
                {"INPUT": self.diskGridPath, "Z_FACTOR": 1.0, "OUTPUT": slope_path},
                feedback=QgsProcessingFeedback(),
            )
        except Exception:
            processing.run(
                "gdal:slope",
                {
                    "INPUT": self.diskGridPath,
                    "BAND": 1,
                    "SCALE": 1.0,
                    "AS_PERCENT": False,
                    "COMPUTE_EDGES": True,
                    "ZEVENBERGEN": False,
                    "OUTPUT": slope_path,
                },
                feedback=QgsProcessingFeedback(),
            )

        slope_layer = QgsRasterLayer(slope_path, self.base_name + "_MRVBF_slope")
        if slope_layer.isValid():
            QgsProject.instance().addMapLayer(slope_layer)

        sds = gdal.Open(slope_path)
        slope_arr = sds.GetRasterBand(1).ReadAsArray().astype(float)
        s_nodata = sds.GetRasterBand(1).GetNoDataValue()
        s_gt = sds.GetGeoTransform()
        sds = None
        if s_nodata is not None:
            slope_arr[slope_arr == s_nodata] = np.nan
        if s_gt[5] > 0:  # south-up -> north-up, matching mrvbf_grid orientation
            slope_arr = np.flipud(slope_arr)
        slope_arr = self._fit_to_shape(slope_arr, mrvbf_grid.shape)

        # --- 3-band composite: R = MRRTF, G = MRVBF, B = slope ---
        rgb_path = self.insert_text_before_extension(self.diskGridPath, "_MRVBF_RGB")
        base, _ext = os.path.splitext(rgb_path)
        rgb_path = base + ".tif"
        if self.is_layer_loaded(self.base_name + "_MRVBF_RGB"):
            project = QgsProject.instance()
            for lyr in project.mapLayersByName(self.base_name + "_MRVBF_RGB"):
                project.removeMapLayer(lyr.id())
        self._write_multiband_raster(
            [mrrtf_grid, mrvbf_grid, slope_arr], rgb_path, self.layer
        )
        rgb_layer = QgsRasterLayer(rgb_path, self.base_name + "_MRVBF_RGB")
        if rgb_layer.isValid():
            QgsProject.instance().addMapLayer(rgb_layer)

        # everything for this operation has already been written and loaded
        self.suffix = ""

    @staticmethod
    def _fit_to_shape(arr, shape):
        """Crop/pad arr (with NaN) so it has exactly `shape`."""
        out = np.full(shape, np.nan, dtype=float)
        ny = min(arr.shape[0], shape[0])
        nx = min(arr.shape[1], shape[1])
        out[:ny, :nx] = arr[:ny, :nx]
        return out

    def _write_multiband_raster(
        self, bands, raster_path, reference_layer, no_data_value=np.nan
    ):
        """Write a list of equal-shaped 2D arrays as a multi-band Float32
        GeoTIFF georeferenced from `reference_layer`."""
        if os.path.exists(raster_path):
            try:
                os.remove(raster_path)
                if os.path.exists(raster_path + ".aux.xml"):
                    os.remove(raster_path + ".aux.xml")
            except OSError:
                self.iface.messageBar().pushMessage(
                    "Couldn't overwrite " + raster_path,
                    level=Qgis.Warning,
                    duration=15,
                )
                return -1

        rows, cols = bands[0].shape
        driver = gdal.GetDriverByName("GTiff")
        output_raster = driver.Create(
            raster_path, cols, rows, len(bands), gdal.GDT_Float32
        )

        provider = reference_layer.dataProvider()
        extent = provider.extent()
        output_raster.SetGeoTransform(
            [
                extent.xMinimum(),
                extent.width() / cols,
                0,
                extent.yMaximum(),
                0,
                -extent.height() / rows,
            ]
        )
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(int(reference_layer.crs().authid().split(":")[1]))
        output_raster.SetProjection(srs.ExportToWkt())

        for i, arr in enumerate(bands, start=1):
            band = output_raster.GetRasterBand(i)
            band.SetNoDataValue(no_data_value)
            band.WriteArray(np.nan_to_num(arr, nan=no_data_value))
            band.FlushCache()
            band = None
        output_raster = None
        return 0

    def procBSDworms(self):
        num_levels = int(self.dlg.spinBox_levels.value())
        bottom_level = self.to_int(self.dlg.doubleSpinBox_base.text())
        delta_z = self.to_float(self.dlg.doubleSpinBox_inc.text())
        layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
        crs = layer.crs()

        if crs.isGeographic():
            self.iface.messageBar().pushMessage(
                "This is a geographic projection, you need to convert it to a projected CRS",
                level=Qgis.Warning,
                duration=15,
            )
            return False
        else:
            if layer.isValid():
                self.diskGridPath = layer.dataProvider().dataSourceUri()
                self.dx = layer.rasterUnitsPerPixelX()
                self.dy = layer.rasterUnitsPerPixelY()
                if not self.validCRS(layer):
                    return False
                crs = int(layer.crs().authid().split(":")[1])

                self.processor = GeophysicalProcessor(self.dx, self.dy, self.buffer)
                shps = self.dlg.checkBox_worms_shp.isChecked()
                if shps:
                    try:
                        import sklearn
                    except ImportError:
                        QMessageBox.information(
                            None,  # Parent widget
                            "",
                            "Missing Packages for SGTool: "  # Window title
                            + f"The following Python packages are required for conversion to shapefile, but not installed: scikit-learn\n\n"
                            "Please open the QGIS Python Console and run the following command:\n\n"
                            f"!pip3 install scikit-learn",  # Message text
                            QMessageBox_Ok,  # Buttons parameter
                        )
                        return False

                """             
                # Access the raster data provider
                provider = layer.dataProvider()

                # Get raster dimensions
                cols = provider.xSize()  # Number of columns
                rows = provider.ySize()  # Number of rows

                # Read raster data as a block
                band = 1  # Specify the band number (1-based index)
                raster_block = provider.block(band, provider.extent(), cols, rows)

                # Copy the block data into a NumPy array
                extent = layer.extent()
                rows, cols = layer.height(), layer.width()
                raster_block = provider.block(1, extent, cols, rows)  # !!!!!
                self.raster_array = np.zeros((rows, cols))
                for i in range(rows):
                    for j in range(cols):
                        self.raster_array[i, j] = raster_block.value(i, j)
            
                """
                # Read raster data via GDAL (consistent row ordering across platforms)
                ds = gdal.Open(self.diskGridPath)
                band = ds.GetRasterBand(1)
                no_data_value = band.GetNoDataValue()  # Band 1
                band_data = ds.GetRasterBand(1).ReadAsArray()  # Always row0=north
                nodata = ds.GetRasterBand(1).GetNoDataValue()
                if nodata is not None:
                    band_data = band_data.astype(float)
                    band_data[band_data == nodata] = np.nan

                self.raster_array = band_data
                if no_data_value is not None:
                    self.raster_array[self.raster_array == no_data_value] = np.nan
                ds = None

                self.processor.bsdwormer(
                    self.raster_array,
                    layer,
                    self.diskGridPath,
                    num_levels,
                    bottom_level,
                    delta_z,
                    shps,
                    crs,
                )
                self.iface.messageBar().pushMessage(
                    "Worms saved to same directory as original grid",
                    level=Qgis.Success,
                    duration=15,
                )

            if shps:

                head_tail = os.path.split(self.diskGridPath)
                out_path = (
                    head_tail[0] + "/" + head_tail[1].split(".")[0] + "_worms.shp"
                )
                layer_name = head_tail[1].split(".")[0] + "_worms.shp"

                layer = QgsVectorLayer(out_path, layer_name)

                if not layer.isValid():
                    raise ValueError(f"Failed to load layer: {out_path}")
                else:
                    # Add the layer to the current QGIS project
                    QgsProject.instance().addMapLayer(layer)

    def set_normalise_in(self):
        """
        Opens a dialog for the user to select an input directory and sets the selected
        directory path to the corresponding line edit widget. If the selected path is
        invalid or empty, displays an error message in the QGIS message bar.
        Functionality:
        - Prompts the user to select a directory using a QFileDialog.
        - Validates the selected directory path.
        - Updates the line edit widget with the valid directory path.
        - Displays an error message if the path is invalid or not selected.
        Raises:
            Displays a critical error message in the QGIS message bar if the path is invalid.
        """

        self.input_directory = QFileDialog.getExistingDirectory(
            None, "Select Input Folder"
        )

        if os.path.exists(self.input_directory) and self.input_directory != "":
            self.dlg.lineEdit_loadPointsPath_normalise_in.setText(self.input_directory)

        else:
            self.iface.messageBar().pushMessage(
                "Error: Path Incorrect",
                level=Qgis.Critical,
                duration=15,
            )

    def set_normalise_out(self):
        """
        Opens a dialog for the user to select an output directory and sets the selected
        directory path to the corresponding line edit widget. If the selected path is invalid
        or empty, an error message is displayed in the QGIS message bar.
        Functionality:
        - Prompts the user to select an output folder using a QFileDialog.
        - Validates the selected directory path.
        - Updates the line edit widget with the selected path if valid.
        - Displays an error message if the path is invalid or empty.
        Raises:
            Displays a critical error message in the QGIS message bar if the selected path
            is invalid or empty.
        """

        self.output_directory = QFileDialog.getExistingDirectory(
            None, "Select Output Folder"
        )

        if os.path.exists(self.output_directory) and self.output_directory != "":
            self.dlg.lineEdit_loadPointsPath_normalise_out.setText(
                self.output_directory
            )

        else:
            self.iface.messageBar().pushMessage(
                "Error: Path Incorrect",
                level=Qgis.Critical,
                duration=15,
            )

    def to_float(self, text):
        """Convert text to float using Qt locale-aware parsing."""
        if isinstance(text, (int, float)):
            return float(text)
        return float(str(text).replace(",", ".").replace(" ", ""))

    def to_int(self, text):
        """Convert text to int using Qt locale-aware parsing."""
        return int(self.to_float(text))

    def util_display_grid(self, grid):
        try:
            import matplotlib.pyplot as plt

        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for some functions, but not installed: matplotlib\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install matplotlib",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False

        plt.imshow(grid, origin="lower", cmap="viridis")
        plt.colorbar(label="Levels")
        plt.title("Grid")
        plt.show()

    def unit_check(self, length):
        """
        Checks if the specified length is valid based on the coordinate reference system (CRS)
        of the selected layer.
        Parameters:
            length (float): The length to be checked.
        Returns:
            bool:
                - False if the CRS is geographic and the length exceeds 100, with a warning
                  message displayed to the user.
                - True otherwise.
        Notes:
            - If the CRS is geographic, lengths should be specified in degrees.
            - A warning message is displayed in the QGIS message bar if the length is invalid.
        """

        selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
        crs = selected_layer.crs()
        if crs.isGeographic() and length > 100:
            self.iface.messageBar().pushMessage(
                "Since this is a geographic projection, you need to specify lengths in degrees",
                level=Qgis.Warning,
                duration=15,
            )
            return False
        else:
            return True

    def _applyAnisoOrientStyle(self):
        layer_name = self.base_name + "_SS_AnisoOrient"
        layers = QgsProject.instance().mapLayersByName(layer_name)
        if not layers:
            return
        layer = layers[0]

        provider = layer.dataProvider()
        stats = provider.bandStatistics(1)

        style = QgsStyle.defaultStyle()
        target = "Full_saturation_spectrum_CCW"
        ramp = style.colorRamp(target)

        if ramp is None:
            # Case-insensitive search in style library
            for name in style.colorRampNames():
                if name.lower() == target.lower():
                    ramp = style.colorRamp(name)
                    break

        if ramp is None:
            # Load from cpt-city catalogue and save into the style library
            cpt_ramp = QgsCptCityColorRamp("ggr/Full_saturation_spectrum_CCW", "")
            if cpt_ramp.count() > 0:
                style.addColorRamp(target, cpt_ramp, True)
                ramp = style.colorRamp(target)
                print(f"Loaded '{target}' from cpt-city and saved to style library.")
            else:
                print(f"Colour ramp '{target}' not found in style library or cpt-city catalogue.")
                return

        ramp_shader = QgsColorRampShader(stats.minimumValue, stats.maximumValue)
        ramp_shader.setColorRampType(QgsColorRampShader.Type.Interpolated)
        ramp_shader.setSourceColorRamp(ramp)
        ramp_shader.classifyColorRamp(classes=256)

        raster_shader = QgsRasterShader()
        raster_shader.setRasterShaderFunction(ramp_shader)

        renderer = QgsSingleBandPseudoColorRenderer(provider, 1, raster_shader)
        layer.setRenderer(renderer)
        layer.setOpacity(0.5)
        layer.triggerRepaint()

    def addNewGrid(self, stdClip=True):
        """
        Adds a new grid layer to the QGIS project. If a layer with the same name already exists,
        it removes the existing layer before adding the new one. The method also handles raster
        file creation, statistics calculation, and renderer configuration for the new layer.
        Steps:
        1. Checks if a layer with the specified name (base_name + suffix) is already loaded.
           If it exists, removes it from the project.
        2. Constructs the file path for the new grid by appending the suffix to the base path.
           Removes any auxiliary XML file associated with the new grid path.
        3. Converts a NumPy array to a raster file and saves it to the constructed path.
        4. Adds the new raster layer to the QGIS project if it is valid.
        5. Calculates statistics for the first band of the raster layer.
        6. Configures the renderer for the raster layer to use contrast enhancement based on
           the calculated statistics.
        Attributes:
            self.suffix (str): Suffix to append to the base name for the new grid layer.
            self.base_name (str): Base name of the grid layer.
            self.diskGridPath (str): File path of the base grid.
            self.new_grid (numpy.ndarray): NumPy array representing the new grid data.
            self.layer (QgsMapLayer): Reference layer for georeferencing the new grid.
            self.diskNewGridPath (str): File path for the new grid with the suffix appended.
        Raises:
            Exception: If the raster file creation fails (indicated by a non-negative error code).
        Notes:
            - The method assumes that the input NumPy array and reference layer are valid.
            - The renderer is configured only if it is of type QgsSingleBandGrayRenderer.
        """

        if self.suffix:
            if self.is_layer_loaded(self.base_name + self.suffix):
                project = QgsProject.instance()
                layer = project.mapLayersByName(self.base_name + self.suffix)[0]
                project.removeMapLayer(layer.id())

            self.diskNewGridPath = self.insert_text_before_extension(
                self.diskGridPath, self.suffix
            )

            if (
                ".tif" not in self.diskNewGridPath.lower()
                and ".tiff" not in self.diskNewGridPath.lower()
            ):
                base, _ = os.path.splitext(self.diskNewGridPath)
                self.diskNewGridPath = base + ".tif"

            if os.path.exists(self.diskNewGridPath + ".aux.xml"):
                os.remove(self.diskNewGridPath + ".aux.xml")
            err = self.numpy_array_to_raster(
                self.new_grid,
                self.diskNewGridPath,
                dx=None,
                xmin=None,
                ymax=None,
                reference_layer=self.layer,
                no_data_value=np.nan,
            )
            if err != -1:
                con_raster_layer = QgsRasterLayer(
                    self.diskNewGridPath, self.base_name + self.suffix
                )
                if con_raster_layer.isValid():
                    # QgsProject.instance().addMapLayer(con_raster_layer)

                    # Add the layer to the project
                    QgsProject.instance().addMapLayer(con_raster_layer)

                    # Access the raster data provider
                    provider = con_raster_layer.dataProvider()

                    # Calculate statistics for the first band
                    band = 1  # Specify the band number
                    stats = provider.bandStatistics(band)

                    # Create or modify the renderer
                    renderer = con_raster_layer.renderer()
                    if isinstance(renderer, QgsSingleBandGrayRenderer):
                        # Set contrast enhancement
                        # Get mean and standard deviation

                        if stdClip:
                            # Calculate min and max values using Mean ± (stddev × 2)
                            mean = stats.mean
                            stddev = stats.stdDev
                            min_value = mean - (stddev * 2)
                            max_value = mean + (stddev * 2)
                            contrast_enhancement = renderer.contrastEnhancement()
                            contrast_enhancement.setMinimumValue(min_value)
                            contrast_enhancement.setMaximumValue(max_value)
                        else:
                            contrast_enhancement = renderer.contrastEnhancement()
                            contrast_enhancement.setMinimumValue(stats.minimumValue)
                            contrast_enhancement.setMaximumValue(stats.maximumValue)

                        # Refresh the layer
                        con_raster_layer.triggerRepaint()
                    else:
                        print("Renderer is not a QgsSingleBandGrayRenderer.")

    def processGeophysics_fft(self):
        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid.currentText()
        self.processGeophysics()

    def processGeophysics_conv(self):
        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid_Conv.currentText()
        self.processGeophysics()

    def processGeophysics(self):
        """
        Processes geophysical data based on the selected raster layer and user-defined parameters.
        This method performs the following steps:
        1. Validates the presence of a raster layer and retrieves its properties.
        2. Reads raster data into a NumPy array, handling NoData values if necessary.
        3. Initializes processing utilities and parameters based on the raster data.
        4. Executes a series of geophysical processing operations based on user selections.
        5. Adds the resulting processed grids to the project.
        Attributes:
            self.localGridName (str): Name of the local raster grid to process.
            self.layer (QgsRasterLayer): The selected raster layer.
            self.base_name (str): Base name of the raster layer.
            self.diskGridPath (str): File path to the raster data source.
            self.dx (float): Raster pixel size in the X direction.
            self.dy (float): Raster pixel size in the Y direction.
            self.raster_array (numpy.ndarray): Array representation of the raster data.
            self.buffer (int): Buffer size for processing, constrained by user input.
        Processing Steps:
            - Directional cleaning
            - Reduction to pole
            - Regional removal
            - Derivative computation
            - Tilt angle calculation
            - Analytic signal computation
            - Upward/downward continuation
            - Band-pass filtering
            - Frequency cut filtering
            - Automatic gain control (AGC)
            - Vertical integration
            - Total horizontal gradient
            - Statistical measures (mean, median, Gaussian filtering)
            - Directional filtering
            - Sun shading
            - Handling NaN values
            - Statistical surface measures (min, max, kurtosis, standard deviation, variance, skewness)
            - Digital terrain model classification
            - Polygon generation
        Post-Processing:
            - Adds new grids to the project after each operation.
            - Resets user interface checkboxes after processing.
        Note:
            This method relies on external classes such as `GeophysicalProcessor`, `ConvolutionFilter`,
            `SG_Util`, and `SpatialStats` for specific processing tasks.
        """

        process = False

        if self.localGridName and self.localGridName != "":
            self.parseParams()

            self.layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
            if self.layer.isValid():
                if not self.validCRS(self.layer):
                    return False
                self.base_name = self.localGridName

                self.diskGridPath = self.layer.dataProvider().dataSourceUri()
                self.dx = abs(self.layer.rasterUnitsPerPixelX())
                self.dy = abs(self.layer.rasterUnitsPerPixelY())

                """
                # Access the raster data provider
                provider = self.layer.dataProvider()

                # Get raster dimensions
                cols = provider.xSize()  # Number of columns
                rows = provider.ySize()  # Number of rows

                # Read raster data as a block
                band = 1  # Specify the band number (1-based index)
                raster_block = provider.block(band, provider.extent(), cols, rows)

                # Copy the block data into a NumPy array
                extent = self.layer.extent()
                rows, cols = self.layer.height(), self.layer.width()
                raster_block = provider.block(1, extent, cols, rows)  # !!!!!
                self.raster_array = np.zeros((rows, cols))
                for i in range(rows):
                    for j in range(cols):
                        self.raster_array[i, j] = raster_block.value(i, j)
                """
                # Read raster data via GDAL (consistent row ordering across platforms)
                ds = gdal.Open(self.diskGridPath)
                band = ds.GetRasterBand(1)
                no_data_value = band.GetNoDataValue()  # call this BEFORE ReadAsArray
                band_data = ds.GetRasterBand(1).ReadAsArray()  # Always row0=north
                if no_data_value is not None:
                    band_data = band_data.astype(float)
                    band_data[band_data == no_data_value] = np.nan
                self.raster_array = band_data
                gt = ds.GetGeoTransform()
                if gt[5] > 0:  # positive pixel height means south-up, flip to north-up
                    self.raster_array = np.flipud(self.raster_array)
                ds = None
                rows, cols = self.raster_array.shape
                process = True

        if process:
            self.buffer = min(rows, cols)
            if self.buffer > self.to_int(self.dlg.lineEdit_13_max_buffer.text()):
                self.buffer = self.to_int(self.dlg.lineEdit_13_max_buffer.text())
            self.processor = GeophysicalProcessor(self.dx, self.dy, self.buffer)

            self.convolution = ConvolutionFilter(self.raster_array)
            self.SG_Util = SG_Util(self.raster_array)
            self.SpatialStats = SpatialStats(self.raster_array)
            self.PCAICA = PCAICA(self.raster_array)

            self.suffix = ""
            if self.DirClean:
                self.procDirClean()
                self.addNewGrid(stdClip=True)
            if self.RTE_P:
                if self.RTE_P_type == "Diff. RTP":
                    self.procVRTP()
                else:
                    self.procRTP_E(sgtool_instance=self)
                self.addNewGrid(stdClip=True)
            if self.RemRegional:
                self.procRemRegional()
                self.addNewGrid(stdClip=True)
            if self.Derivative:
                self.procDerivative()
                self.addNewGrid(stdClip=True)
            if self.TA:
                self.procTiltAngle()
                self.addNewGrid(stdClip=True)
            if self.AS:
                self.procAnalyticSignal()
                self.addNewGrid(stdClip=True)
            if self.Continuation:
                self.procContinuation()
                self.addNewGrid(stdClip=True)
            if self.BandPass:
                self.procBandPass()
                self.addNewGrid(stdClip=True)
            if self.FreqCut:
                self.procFreqCut()
                self.addNewGrid(stdClip=True)
            if self.AGC:
                self.procAGC()
                self.addNewGrid(stdClip=True)
            if self.VI:
                self.procvInt()
                self.addNewGrid(stdClip=True)
            if self.THG:
                self.procTHG()
                self.addNewGrid(stdClip=True)

            if self.Mean:
                self.procMean()
                self.addNewGrid(stdClip=True)
            if self.Median:
                self.procMedian()
                self.addNewGrid(stdClip=True)
            if self.Gaussian:
                self.procGaussian()
                self.addNewGrid(stdClip=True)
            if self.Direction:
                self.procDirectional()
                self.addNewGrid(stdClip=True)
            if self.SunShade:
                self.procSunShade()
                self.addNewGrid(stdClip=False)
            if self.NaN:
                self.procNaN()
                self.addNewGrid(stdClip=True)

            if self.SS_Min:
                self.procSS_Min()
                self.addNewGrid(stdClip=True)
            if self.SS_Max:
                self.procSS_Max()
                self.addNewGrid(stdClip=True)
            if self.SS_Kurtosis:
                self.procSS_Kurtosis()
                self.addNewGrid(stdClip=True)
            if self.SS_StdDev:
                self.procSS_StdDev()
                self.addNewGrid(stdClip=True)
            if self.SS_Variance:
                self.procSS_Variance()
                self.addNewGrid(stdClip=True)
            if self.SS_Skewness:
                self.procSS_Skewness()
                self.addNewGrid(stdClip=True)
            if self.SS_Anisotropy:
                self.procSS_Anisotropy()
                self.addNewGrid(stdClip=True)   # adds _SS_AnisoMag
                self.new_grid = self._aniso_orientation
                self.suffix = "_SS_AnisoOrient"
                self.addNewGrid(stdClip=False)  # adds _SS_AnisoOrient
                self._applyAnisoOrientStyle()
            if self.SS_ChainLength:
                self.procSS_ChainLength()
            if self.SS_Streamline:
                self.procSS_Streamline()
            if self.MRVBF:
                self.procMRVBF()
            if self.PCA:
                self.procPCA()
            if self.ICA:
                self.procICA()
            if self.Polygons:
                self.procPolygons()
            if self.ED:
                self.procEulerDeconvolution(self.raster_array)

            self.resetCheckBoxes()

    def resetCheckBoxes(self):
        """
        Resets all checkboxes and associated boolean flags in the dialog to their default unchecked/False state.
        This method performs the following actions:
        - Sets all checkboxes in the dialog to an unchecked state using `setChecked(False)`.
        - Resets all associated boolean flags to `False`.
        Checkboxes reset include:
        - Processing options (e.g., RTE_P, tilt derivative, analytic signal, etc.).
        - Filtering options (e.g., mean, median, Gaussian, directional, etc.).
        - Statistical options (e.g., minimum, maximum, kurtosis, etc.).
        - Miscellaneous options (e.g., NaN handling, polygons, etc.).
        Boolean flags reset include:
        - Processing flags (e.g., RTE_P, TA, AS, etc.).
        - Filtering flags (e.g., Mean, Median, Gaussian, etc.).
        - Miscellaneous flags (e.g., Polygons, SunShade, etc.).
        """

        self.dlg.checkBox_4_RTE_P.setChecked(False)
        self.dlg.checkBox_7_tiltDerivative.setChecked(False)
        self.dlg.checkBox_8_analyticSignal.setChecked(False)
        self.dlg.checkBox_4_PGrav.setChecked(False)
        self.dlg.checkBox_9_continuation.setChecked(False)

        self.dlg.checkBox_3_DirClean.setChecked(False)
        self.dlg.checkBox_5_regional.setChecked(False)
        self.dlg.checkBox_10_bandPass.setChecked(False)
        self.dlg.checkBox_10_freqCut.setChecked(False)
        self.dlg.checkBox_11_1vd_agc.setChecked(False)

        self.dlg.checkBox_6_derivative.setChecked(False)
        self.dlg.checkBox_11_tot_hz_grad.setChecked(False)

        self.dlg.checkBox_Mean.setChecked(False)
        self.dlg.checkBox_Median.setChecked(False)
        self.dlg.checkBox_Gaussian.setChecked(False)
        self.dlg.checkBox_Directional.setChecked(False)
        self.dlg.checkBox_SunShading.setChecked(False)

        self.dlg.checkBox_NaN.setChecked(False)
        self.dlg.checkBox_polygons.setChecked(False)

        self.dlg.checkBox_SS_Min.setChecked(False)
        self.dlg.checkBox_SS_Max.setChecked(False)
        self.dlg.checkBox_SS_Kurtosis.setChecked(False)
        self.dlg.checkBox_SS_StdDev.setChecked(False)
        self.dlg.checkBox_SS_Variance.setChecked(False)
        self.dlg.checkBox_SS_Skewness.setChecked(False)
        self.dlg.checkBox_SS_Anisotropy.setChecked(False)
        self.dlg.checkBox_SS_ChainLength.setChecked(False)
        self.dlg.checkBox_SS_Streamline.setChecked(False)
        self.dlg.checkBox_MRVBF.setChecked(False)
        self.dlg.checkBox_PCA.setChecked(False)
        self.dlg.checkBox_ICA.setChecked(False)
        self.dlg.checkBox_ED_Stats.setChecked(False)
        self.dlg.checkBox_ED.setChecked(False)

        self.RTE_P = False
        self.TA = False
        self.AS = False

        self.DirClean = False
        self.RemRegional = False
        self.Derivative = False
        self.Continuation = False
        self.BandPass = False
        self.AGC = False
        self.FreqCut = False
        self.VI = False
        self.THG = False
        self.Mean = False
        self.Median = False
        self.Gaussian = False
        self.Gaussian = False
        self.Direction = False
        self.SunShade = False
        self.Polygons = False

    def is_layer_loaded(self, layer_name):
        """
        Check if a layer with the specified name is already loaded in QGIS.

        Parameters:
            layer_name (str): The name of the layer to check.

        Returns:
            bool: True if the layer is loaded, False otherwise.
        """
        for layer in QgsProject.instance().mapLayers().values():
            if layer.name() == layer_name:
                return True
        return False

    def select_RGBgrid_file(self):
        start_directory = self.last_directory if self.last_directory else os.getcwd()

        self.diskRGBGridPath, _filter = QFileDialog.getOpenFileName(
            None,
            "Select RGB Image File",
            start_directory,
            "Grids (*.TIF *.tif *.TIFF *.tiff)",
        )
        if os.path.exists(self.diskRGBGridPath) and self.diskRGBGridPath != "":
            self.dlg.lineEdit_2_loadGridPath_2.setText(self.diskRGBGridPath)
            self.last_directory = os.path.dirname(self.diskRGBGridPath)
            basename = os.path.basename(self.diskRGBGridPath)
            filename_without_extension = os.path.splitext(basename)[0]
            self.rgb_layer = QgsRasterLayer(
                self.diskRGBGridPath, filename_without_extension
            )
            if not self.is_layer_loaded(self.diskRGBGridPath):
                QgsProject.instance().addMapLayer(self.rgb_layer)

    def _on_pick_rgb_clicked(self):
        canvas = self.iface.mapCanvas()
        self._prev_map_tool = canvas.mapTool()
        self._rgb_picker = RGBPickerMapTool(canvas)
        self._rgb_picker.canvasClicked.connect(self._on_rgb_picked)
        canvas.setMapTool(self._rgb_picker)

    def _on_rgb_picked(self, point, button):
        from qgis.PyQt.QtCore import QTimer
        canvas = self.iface.mapCanvas()
        prev, picker = self._prev_map_tool, self._rgb_picker

        def _restore():
            if prev:
                canvas.setMapTool(prev)
            else:
                canvas.unsetMapTool(picker)

        # Defer tool restoration until after the current event cycle so QGIS
        # cannot re-dispatch this mouse event to the restored (identify) tool.
        QTimer.singleShot(0, _restore)

        rgb_path = self.dlg.lineEdit_2_loadGridPath_2.text()
        if not rgb_path or not os.path.exists(rgb_path):
            self.iface.messageBar().pushMessage(
                "SGTool", "No RGB grid loaded", level=Qgis.Warning, duration=5
            )
            return

        ds = gdal.Open(rgb_path)
        if ds is None or ds.RasterCount < 3:
            self.iface.messageBar().pushMessage(
                "SGTool", "Grid must have at least 3 bands", level=Qgis.Warning, duration=5
            )
            return

        # Transform canvas point to layer CRS if needed
        layer_crs = QgsCoordinateReferenceSystem()
        layer_crs.createFromWkt(ds.GetProjection())
        canvas_crs = canvas.mapSettings().destinationCrs()
        sample_point = point
        if layer_crs.isValid() and layer_crs != canvas_crs:
            try:
                transform = QgsCoordinateTransform(
                    canvas_crs, layer_crs, QgsProject.instance()
                )
                sample_point = transform.transform(point)
            except Exception as e:
                print(f"SGTool: CRS transform failed, using untransformed point: {e}")

        inv_gt = gdal.InvGeoTransform(ds.GetGeoTransform())
        px, py = gdal.ApplyGeoTransform(inv_gt, sample_point.x(), sample_point.y())
        px, py = int(px), int(py)

        if px < 0 or py < 0 or px >= ds.RasterXSize or py >= ds.RasterYSize:
            self.iface.messageBar().pushMessage(
                "SGTool", "Click is outside image extent", level=Qgis.Warning, duration=5
            )
            return

        r = int(ds.GetRasterBand(1).ReadAsArray(px, py, 1, 1)[0][0])
        g = int(ds.GetRasterBand(2).ReadAsArray(px, py, 1, 1)[0][0])
        b = int(ds.GetRasterBand(3).ReadAsArray(px, py, 1, 1)[0][0])

        current = self.dlg.textEdit_2_colour_list.toPlainText().strip()
        new_triplet = f"{r}, {g}, {b}"
        if current:
            new_text = current + ", " + new_triplet
        else:
            new_text = new_triplet
        self.dlg.textEdit_2_colour_list.setPlainText(new_text)

    def processRGB(self):
        """
        Processes an RGB grid file by converting it to grayscale using a specified CSS color list
        and loads the resulting raster layer into the QGIS project.
        This method performs the following steps:
        1. Reads the file path of the RGB grid from the user interface.
        2. Validates the existence of the file and the presence of a CSS color list.
        3. Converts the RGB grid to grayscale using the provided CSS color list.
        4. Loads the resulting grayscale raster layer into the QGIS project if it is not already loaded.
        Raises appropriate messages in the QGIS message bar for the following cases:
        - If the file path is empty or does not exist.
        - If the CSS color list is not defined.
        - If the conversion fails due to invalid CSS color names.
        Attributes:
            self.diskRGBGridPath (str): The file path of the RGB grid.
            self.layer (QgsRasterLayer): The resulting grayscale raster layer.
        Returns:
            None
        """

        self.diskRGBGridPath = self.dlg.lineEdit_2_loadGridPath_2.text()
        if self.diskRGBGridPath != "":
            if os.path.exists(self.diskRGBGridPath):
                LUT_list = self.dlg.textEdit_2_colour_list.toPlainText()

                if LUT_list != "":

                    result, RGBGridPath_gray = self.convert_RGB_to_grey(
                        self.diskRGBGridPath, LUT_list
                    )
                    if result:

                        basename = os.path.basename(RGBGridPath_gray)
                        filename_without_extension = os.path.splitext(basename)[0]

                        self.layer = QgsRasterLayer(
                            RGBGridPath_gray, filename_without_extension
                        )
                        """try:
                            test_proj = self.layer.crs().authid()
                            self.layer.setCrs(test_proj)

                        except:
                            # Define the new CRS (e.g., EPSG:4326 for WGS84)
                            new_crs = QgsCoordinateReferenceSystem("EPSG:4326")
                            # Set the CRS for the raster layer
                            self.layer.setCrs(new_crs)"""
                        if not self.is_layer_loaded(filename_without_extension):
                            QgsProject.instance().addMapLayer(self.layer)

                    else:
                        if RGBGridPath_gray != -3:
                            self.iface.messageBar().pushMessage(
                                "Conversion failed, check CSS colour names",
                                level=Qgis.Warning,
                                duration=15,
                            )

                else:
                    self.iface.messageBar().pushMessage(
                        "First define a CSS Colour list <a href='https://matplotlib.org/stable/gallery/color/named_colors.html#css-colors'> (See here for list of colours)</a>",
                        level=Qgis.Info,
                        duration=15,
                    )

    def select_grid_file(self):
        """
        Opens a file dialog to select a grid file, processes the selected file, and updates the UI accordingly.
        This method allows the user to select a grid file from the file system. It supports multiple file formats
        such as TIF, TIFF, GRD, and ERS. Depending on the file type, it performs specific operations like extracting
        the CRS (Coordinate Reference System) or loading the file as a raster layer in QGIS.
        Steps:
        1. Opens a file dialog starting from the last used directory or the current working directory.
        2. Updates the UI with the selected file path and enables processing options.
        3. Handles specific file types:
           - For `.grd` files, attempts to extract the CRS from an accompanying XML file. Defaults to EPSG:4326 if no CRS is found.
           - For `.tif` and `.ers` files, loads the file as a raster layer in QGIS.
        UI Updates:
        - Updates the text field with the selected file path.
        - Enables the "Apply Processing" button if a valid file is selected.
        - Displays messages in the QGIS message bar for CRS-related information.
        Raises:
        - Displays warnings in the QGIS message bar if no CRS is found for `.grd` files.
        Attributes:
        - self.diskGridPath: Stores the path of the selected grid file.
        - self.last_directory: Updates to the directory of the selected file for future use.
        - self.layer: Stores the loaded raster layer for `.tif` and `.ers` files.
        File Filters:
        - Supports the following file extensions (upper or lower case) :  `.tif`,  `.tiff`, `.grd`,  `.ers`, 'mag', 'grv'
        Note:
        - The method assumes the presence of an XML file for `.grd` files to extract CRS information.
        - If the layer is not already loaded, it adds the raster layer to the QGIS project.
        """

        start_directory = self.last_directory if self.last_directory else os.getcwd()
        file_filter = "Grids (*.TIF *.tif *.TIFF *.tiff *.grd *.GRD *.ERS *.ers *.mag *.MAG *.grv *.GRV)"

        self.diskGridPath, _filter = QFileDialog.getOpenFileName(
            None,
            "Select Data File",
            start_directory,
            file_filter,  # "Grids (*.TIF;*.tif;*.TIFF;*.tiff;*.grd;*GRD;*.ERS;*.ers)",
        )
        suffix = self.diskGridPath.split(".")[-1].lower()
        epsg = None
        if os.path.exists(self.diskGridPath) and self.diskGridPath != "":
            self.dlg.lineEdit_2_loadGridPath.setText(self.diskGridPath)
            self.dlg.pushButton_3_applyProcessing.setEnabled(True)
            self.last_directory = os.path.dirname(self.diskGridPath)
            if suffix == "grd":
                if os.path.exists(self.diskGridPath + ".xml"):
                    epsg = extract_proj_str(self.diskGridPath + ".xml")
                if epsg == None:
                    epsg = 4326
                    self.iface.messageBar().pushMessage(
                        "No CRS found in XML, default to 4326",
                        level=Qgis.Warning,
                        duration=15,
                    )
                else:
                    self.iface.messageBar().pushMessage(
                        "CRS Read from XML as " + str(epsg),
                        level=Qgis.Info,
                        duration=15,
                    )
                # self.dlg.mQgsProjectionSelectionWidget.setCrs(QgsCoordinateReferenceSystem('EPSG:'+str(epsg)))
                self.save_a_grid(epsg)
            elif suffix in ("tif", "tiff", "ers"):
                basename = os.path.basename(self.diskGridPath)
                filename_without_extension = os.path.splitext(basename)[0]
                self.layer = QgsRasterLayer(
                    self.diskGridPath, filename_without_extension
                )

                if not self.is_layer_loaded(self.diskGridPath):
                    QgsProject.instance().addMapLayer(self.layer)
            elif suffix == "mag" or suffix == "grv":
                self.loadNoddyGrid(self.diskGridPath)

    def loadNoddyGrid(self, input_file=None):
        """
        Reads an ASCII file with magnetic/gravity data and converts it to GeoTIFF format.
        Also creates a text file with inclination and declination information.

        Parameters:
        input_file (str): Path to the input ASCII file

        Returns:
        tuple: (geotiff_path, ntxt_path) - paths to the created files
        """

        # Read the file
        with open(input_file, "r") as f:
            lines = f.readlines()

        # Find the start of the header data by skipping comment lines
        header_start = 0
        for i, line in enumerate(lines):
            if not line.strip().startswith("#") and line.strip():
                header_start = i
                break

        # Parse header information (6 lines after file code and comments)
        file_code = int(lines[header_start].strip())
        if file_code == 333:
            dataType = "mag"
        else:
            dataType = "grv"

        # Skip any additional comment lines after file code
        data_line_index = header_start + 1
        while data_line_index < len(lines) and lines[
            data_line_index
        ].strip().startswith("#"):
            data_line_index += 1

        # Now read the 6 header lines
        calc_range, cols, rows, layers = map(
            int, lines[data_line_index].strip().split()
        )
        inclination, declination, intensity = map(
            float, lines[data_line_index + 1].strip().split()
        )
        sw_x, sw_y, sw_z = map(float, lines[data_line_index + 2].strip().split())
        ne_x, ne_y, ne_z = map(float, lines[data_line_index + 3].strip().split())
        cell_size, survey_height = map(
            float, lines[data_line_index + 4].strip().split()
        )

        # load grid data
        data_array = mag = np.loadtxt(input_file, skiprows=8)
        actual_rows = mag.shape[0]
        actual_cols = mag.shape[1]

        # Calculate geotransform parameters
        # [top_left_x, pixel_width, rotation, top_left_y, rotation, pixel_height]
        pixel_width = (ne_x - sw_x) / actual_cols
        pixel_height = (
            -(ne_y - sw_y) / actual_rows
        )  # Negative because Y decreases from north to south

        # Top-left corner coordinates (GDAL convention)
        top_left_x = sw_x
        top_left_y = ne_y  # Start from the north (top)

        geotransform = [top_left_x, pixel_width, 0, top_left_y, 0, pixel_height]

        # Set up file paths
        base_name = os.path.splitext(input_file)[0]
        geotiff_path = f"{base_name}_{dataType}.tif"
        head_tail = os.path.split(input_file)
        layer_name = head_tail[1].split(".")[0] + f"_{dataType}"

        if os.path.exists(geotiff_path):
            try:
                os.remove(geotiff_path)
                if os.path.exists(geotiff_path + "aux.xml"):
                    os.remove(geotiff_path + "aux.xml")
            except:
                self.iface.messageBar().pushMessage(
                    "Couldn't delete layer, may be open in another program? On windows files on non-C: drive may be hard to delete",
                    level=Qgis.Warning,
                    duration=15,
                )
                return -1

        # Create the GeoTIFF using GDAL
        driver = gdal.GetDriverByName("GTiff")

        # Create the dataset
        dataset = driver.Create(
            geotiff_path,
            actual_cols,
            actual_rows,
            1,
            gdal.GDT_Float32,  # Number of bands
        )

        if dataset is None:
            raise RuntimeError(f"Could not create {geotiff_path}")

        # Set the geotransform
        dataset.SetGeoTransform(geotransform)

        # Set the projection (UTM29N)
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(32629)  # UTM29N
        dataset.SetProjection(srs.ExportToWkt())

        # Write the data to the band
        band = dataset.GetRasterBand(1)
        band.WriteArray(data_array)

        # Set nodata value if needed
        # band.SetNoDataValue(-9999)

        # Set metadata
        metadata = {
            "file_code": str(file_code),
            "calculation_range": str(calc_range),
            "original_rows": str(rows),
            "original_cols": str(cols),
            "layers": str(layers),
            "intensity": str(intensity),
            "inclination": str(inclination),
            "declination": str(declination),
            "sw_corner_x": str(sw_x),
            "sw_corner_y": str(sw_y),
            "sw_corner_z": str(sw_z),
            "ne_corner_x": str(ne_x),
            "ne_corner_y": str(ne_y),
            "ne_corner_z": str(ne_z),
            "cube_size": str(cell_size),
            "survey_height": str(survey_height),
        }
        dataset.SetMetadata(metadata)

        # Properly close the dataset
        dataset.FlushCache()
        dataset = None

        self.layer = QgsRasterLayer(geotiff_path, layer_name)
        if not self.is_layer_loaded(layer_name):
            QgsProject.instance().addMapLayer(self.layer)

    # save grd file as geotiff
    def save_a_grid(self, epsg):
        """
        Saves a grid file as a GeoTIFF and loads it as a raster layer in QGIS.

        This method processes a grid file, converts it to a GeoTIFF format, and
        optionally loads it into the QGIS project. It handles various file
        operations, including checking for file existence, removing old files,
        and setting geospatial metadata.

        Args:
            epsg (int): The EPSG code representing the spatial reference system
                        to be used for the output GeoTIFF.

        Raises:
            None. However, it displays warning messages in the QGIS message bar
            if issues are encountered, such as missing files or unsupported data types.

        Notes:
            - The method assumes that `self.diskGridPath` contains the path to the
              input grid file.
            - If the input grid file is not found or not selected, a warning is
              displayed, and the method exits.
            - The method uses GDAL for file creation and geospatial metadata
              assignment.
            - The output GeoTIFF file is saved in the same directory as the input
              grid file, with the same base name but a `.tif` extension.
            - If the output file already exists, it is removed before creating a
              new one.
            - The method supports only certain data types and will display a
              warning if unsupported types are encountered.
            - The resulting raster layer is added to the QGIS project if it is
              not already loaded.
        """

        # load grd file and store in memory
        if self.diskGridPath != "":
            if not os.path.exists(self.diskGridPath):
                self.iface.messageBar().pushMessage(
                    "File: " + self.diskGridPath + " not found",
                    level=Qgis.Warning,
                    duration=3,
                )
            else:
                grid, header, Gdata_type = load_oasis_montaj_grid_optimized(
                    self.diskGridPath
                )
                # grid,header,Gdata_type=load_oasis_montaj_grid(self.diskGridPath)
                if Gdata_type == -1:
                    self.iface.messageBar().pushMessage(
                        "Sorry, can't read 'SHORT' or 'INT' data types at the moment",
                        level=Qgis.Warning,
                        duration=30,
                    )
                    return
                else:
                    directory_path = os.path.dirname(self.diskGridPath)
                    basename = os.path.basename(self.diskGridPath)
                    filename_without_extension = os.path.splitext(basename)[0]
                    self.diskGridPath = (
                        directory_path + "/" + filename_without_extension + ".tif"
                    )

                    fn = self.diskGridPath
                    if os.path.exists(self.diskGridPath) and not self.is_layer_loaded(
                        filename_without_extension
                    ):
                        os.remove(self.diskGridPath)
                        if os.path.exists(self.diskGridPath + ".aux.xml"):
                            os.remove(self.diskGridPath + ".aux.xml")

                    basename = os.path.basename(self.diskGridPath)
                    extension = os.path.splitext(basename)[1].lower()
                    if extension == "ers":
                        driver = gdal.GetDriverByName("ERS")
                    else:
                        driver = gdal.GetDriverByName("GTiff")

                    if header["ordering"] == 1:
                        ds = driver.Create(
                            fn,
                            xsize=header["shape_e"],
                            ysize=header["shape_v"],
                            bands=1,
                            eType=Gdata_type,
                        )
                    else:
                        ds = driver.Create(
                            fn,
                            xsize=header["shape_v"],
                            ysize=header["shape_e"],
                            bands=1,
                            eType=Gdata_type,
                        )

                    ds.GetRasterBand(1).WriteArray(grid)
                    geot = [
                        header["x_origin"] - (header["spacing_e"] / 2),
                        header["spacing_e"],
                        0,
                        header["y_origin"] - (header["spacing_v"] / 2),
                        0,
                        header["spacing_e"],
                    ]
                    ds.SetGeoTransform(geot)
                    srs = osr.SpatialReference()

                    srs.ImportFromEPSG(int(epsg))
                    ds.SetProjection(srs.ExportToWkt())
                    ds.FlushCache()
                    ds = None

                    self.layer = QgsRasterLayer(
                        self.diskGridPath, filename_without_extension
                    )
                    if not self.is_layer_loaded(filename_without_extension):
                        QgsProject.instance().addMapLayer(self.layer)

        else:
            self.iface.messageBar().pushMessage(
                "You need to select a file first", level=Qgis.Warning, duration=3
            )

    def select_point_file(self):
        start_directory = self.last_directory if self.last_directory else os.getcwd()
        file_filter = (
            "points or lines (*.csv *.txt *.dat *.xyz *.CSV *.DAT *.TXT *.XYZ)"
        )

        self.diskPointsPath, _filter = QFileDialog.getOpenFileName(
            None,
            "Select Data File",
            start_directory,
            file_filter,  # "points or lines (*.csv;*.txt;*.xyz;*.CSV;*.TXT;*.XYZ)",
        )
        if os.path.exists(self.diskPointsPath) and self.diskPointsPath != "":
            self.last_directory = os.path.dirname(self.diskPointsPath)

            basename = os.path.basename(self.diskPointsPath)
            extension = os.path.splitext(basename)[1]
            self.dlg.lineEdit_loadPointsPath.setText(self.diskPointsPath)

            if extension.upper() == ".XYZ":
                self.line_data_cols = self.get_XYZ_header(self.diskPointsPath)
                self.dlg.mQgsProjectionSelectionWidget.setEnabled(True)
                self.dlg.pushButton_load_point_data.setEnabled(True)
                self.pointType = "line"

            elif extension.upper() == ".TXT" or extension.upper() == ".CSV":
                columns = self.read_header(self.diskPointsPath, ",")
                # columns = list(points.columns)
                self.dlg.comboBox_grid_x.setEnabled(True)
                self.dlg.comboBox_grid_y.setEnabled(True)
                self.dlg.mQgsProjectionSelectionWidget.setEnabled(True)
                self.dlg.comboBox_grid_x.clear()
                self.dlg.comboBox_grid_y.clear()

                self.dlg.comboBox_grid_x.addItems(columns)
                self.dlg.comboBox_grid_y.addItems(columns)
                self.dlg.pushButton_load_point_data.setEnabled(True)
                self.pointType = "point"
            else:  # DAT format
                self.parser = AsegGdf2Parser()

                directory_path = os.path.dirname(self.diskPointsPath)
                basename = os.path.basename(self.diskPointsPath)
                filename_without_extension = os.path.splitext(basename)[0]
                extension = os.path.splitext(basename)[1]
                if extension == ".dat":
                    dfnPath = directory_path + "/" + filename_without_extension + ".dfn"
                else:
                    dfnPath = directory_path + "/" + filename_without_extension + ".DFN"

                self.header_list, self.field_defs, self.points_epsg = (
                    self.parser.parse_dfn_file(dfnPath)
                )
                print("self.header_list", self.header_list)
                if self.header_list is None:
                    return

                self.dlg.comboBox_grid_x.setEnabled(True)
                self.dlg.comboBox_grid_y.setEnabled(True)
                self.dlg.mQgsProjectionSelectionWidget.setEnabled(True)

                # Create a CRS object with your specific EPSG code

                if self.points_epsg is None:
                    crs = QgsCoordinateReferenceSystem("EPSG:4326")
                else:
                    crs = QgsCoordinateReferenceSystem(f"EPSG:{self.points_epsg}")

                # Set the CRS on the projection selection widget
                self.dlg.mQgsProjectionSelectionWidget.setCrs(crs)
                self.dlg.comboBox_grid_x.clear()
                self.dlg.comboBox_grid_y.clear()
                self.dlg.comboBox_grid_x.addItems(self.header_list)
                self.dlg.comboBox_grid_y.addItems(self.header_list)
                self.dlg.pushButton_load_point_data.setEnabled(True)
                self.pointType = "point"

    def get_XYZ_header(self, csv_file):
        # Input file path
        # csv_file = r"//wsl.localhost/Ubuntu-20.04/home/mark/gridding/MAG.XYZ"  # Replace with your actual file path

        # Initialize variables
        data_list = []
        current_line_number = None

        # Read the file line-by-line
        with open(csv_file, "r") as file:
            for line in file:
                line = line.strip()
                if line.startswith("LINE:"):  # Check for 'LINE:' markers
                    current_line_number = int(
                        re.search(r"\d+", line).group()
                    )  # Extract the line number
                elif current_line_number is not None:
                    try:
                        # Parse numerical lines
                        parts = list(map(float, line.split()))
                        if len(parts) >= 3:  # Ensure at least 5 components (x, y, z)
                            data_list.append(parts + [current_line_number])
                            return len(parts) - 2
                    except ValueError:
                        pass

    def read_header(self, file_path, delimiter=","):
        """
        Reads the first line of a CSV file and stores the values as a list.

        :param file_path: Path to the CSV file
        :return: List of header values
        """
        with open(file_path, "r") as file:
            header = file.readline().strip().split(delimiter)
            header = list(filter(None, header))
        return header

    def numpy_array_to_raster(
        self,
        numpy_array,
        raster_path,
        dx=None,
        xmin=None,
        ymax=None,
        reference_layer=None,
        no_data_value=np.nan,
    ):
        """
        Convert a NumPy array to a GeoTIFF raster file.

        Parameters:
            numpy_array (numpy.ndarray): The NumPy array to convert.
            raster_path (str): The path to save the raster file.
            reference_layer (QgsRasterLayer, optional): A reference layer for CRS and geotransform.
            no_data_value: Value to use for no data (default is NaN).
        """

        # Check if the file already exists and remove it
        if os.path.exists(raster_path):
            try:
                os.remove(raster_path)
                if os.path.exists(raster_path + "aux.xml"):
                    os.remove(raster_path + "aux.xml")
            except:
                self.iface.messageBar().pushMessage(
                    "Couldn't delete layer, may be open in another program? On windows files on non-C: drive may be hard to delete",
                    level=Qgis.Warning,
                    duration=15,
                )
                return -1

        rows, cols = numpy_array.shape
        driver = gdal.GetDriverByName("GTiff")
        output_raster = driver.Create(raster_path, cols, rows, 1, gdal.GDT_Float32)

        # Set geotransform and projection if a reference layer is provided
        if reference_layer:
            provider = reference_layer.dataProvider()
            extent = provider.extent()
            geotransform = [
                extent.xMinimum(),
                extent.width() / cols,  # pixel width
                0,
                extent.yMaximum(),
                0,
                -extent.height() / rows,  # pixel height (negative)
            ]
            output_raster.SetGeoTransform(geotransform)

            # Set CRS
            srs = osr.SpatialReference()
            epsg_code = (
                reference_layer.crs().authid().split(":")[1]
            )  # Extract the EPSG number
            srs.ImportFromEPSG(int(epsg_code))
            output_raster.SetProjection(srs.ExportToWkt())
        else:
            crs = self.dlg.mQgsProjectionSelectionWidget.crs().authid()
            srs = osr.SpatialReference()
            srs.ImportFromEPSG(int(crs.split(":")[1]))
            output_raster.SetProjection(srs.ExportToWkt())
            geotransform = [
                xmin,
                dx,  # pixel width
                0,
                ymax,
                0,
                -abs(-dx),  # pixel height (negative)
            ]
            output_raster.SetGeoTransform(geotransform)

        # Write data to raster
        band = output_raster.GetRasterBand(1)
        if no_data_value is not None:
            band.SetNoDataValue(no_data_value)
        numpy_array = np.nan_to_num(
            numpy_array, nan=no_data_value
        )  # Replace NaN with no_data_value
        band.WriteArray(numpy_array)
        band.FlushCache()
        # Release band before closing dataset — on Windows, keeping a band
        # reference alive prevents the GDAL dataset from being truly closed,
        # which causes a file-lock conflict when QGIS immediately tries to
        # reopen the same file as a raster layer (most visible when two
        # addNewGrid calls happen back-to-back with no computation in between).
        band = None
        output_raster = None  # Close the file
        return 0

    def day_month_to_decimal_year(self, year, month, day):
        """
        Convert a day and month into a decimal year.

        Parameters:
            year (int): The year.
            month (int): The month (1-12).
            day (int): The day (1-31 depending on the month).

        Returns:
            float: The decimal year.
        """
        # Create datetime object for the given date
        date = datetime(year, month, day)

        # Calculate the start and end of the year
        start_of_year = datetime(year, 1, 1)
        end_of_year = datetime(year + 1, 1, 1)

        # Calculate the total days in the year
        days_in_year = (end_of_year - start_of_year).days

        # Calculate the number of days since the start of the year
        days_since_start_of_year = (date - start_of_year).days

        # Compute the decimal year
        decimal_year = year + days_since_start_of_year / days_in_year

        return decimal_year

    def get_grid_centroid(self, layer):
        """
        Calculate the centroid (midpoint) of the grid based on the extent of the given raster layer.
        Args:
            layer (QgsRasterLayer): The raster layer for which the grid centroid is to be calculated.
        Returns:
            tuple: A tuple containing the x and y coordinates of the grid centroid (midx, midy).
        Note:
            The method assumes that the input layer is valid. If the layer is invalid,
            the behavior of the method is undefined.
        """

        if layer.isValid():
            extent = layer.extent()  # Get the extent of the raster layer

            # calculate midpoint of grid
            midx = extent.xMinimum() + (extent.xMaximum() - extent.xMinimum()) / 2
            midy = extent.yMinimum() + (extent.yMaximum() - extent.yMinimum()) / 2

        return midx, midy

    def validCRS(self, layer):
        """
        Check if the given layer has a valid Coordinate Reference System (CRS).

        Parameters:
            layer (QgsRasterLayer or QgsVectorLayer): The layer to check.

        Returns:
            bool: True if the layer has a valid CRS, False otherwise.
        """
        if layer.crs().authid():
            test_crs = layer.crs().authid().split(":")
            if len(test_crs) != 2:
                valid = False
            else:
                valid = True
        else:
            valid = False

        if not valid:
            self.iface.messageBar().pushMessage(
                "This layer has an unrecognised coordinate reference system (CRS). Please set a valid CRS.",
                level=Qgis.Warning,
                duration=30,
            )
        return valid

    # estimate mag field from centroid of data, date and sensor height
    def update_mag_field(self):
        """
        Updates the magnetic field parameters based on the selected grid layer and user inputs.

        This method retrieves the selected grid layer, calculates the magnetic field parameters
        (inclination, declination, and intensity) using the IGRF model, and updates the corresponding
        widgets in the user interface.

        Steps:
        1. Loads the selected grid layer if it exists.
        2. Retrieves user inputs for magnetic intensity and survey date.
        3. Calculates the centroid of the grid layer and converts it to latitude/longitude.
        4. Computes the magnetic field parameters using the IGRF model.
        5. Updates the user interface with the calculated values.

        Raises a warning if the projection system of the grid layer cannot be interpreted.

        Attributes:
            self.localGridName (str): Name of the selected grid layer.
            self.diskGridPath (str): Path to the grid file on disk.
            self.layer (QgsVectorLayer): The selected grid layer.
            self.base_name (str): Base name of the grid layer.
            self.magn_int (str): Magnetic intensity input by the user.
            self.magn_SurveyDay (int): Day of the survey date.
            self.magn_SurveyMonth (int): Month of the survey date.
            self.magn_SurveyYear (int): Year of the survey date.
            self.RTE_P_inc (float): Calculated magnetic inclination.
            self.RTE_P_dec (float): Calculated magnetic declination.
            self.RTE_P_int (float): Calculated magnetic intensity.

        Notes:
            - Requires the `pyproj` library for coordinate transformations.
            - Assumes the grid layer has a valid CRS (Coordinate Reference System).
            - Displays a warning message if the CRS is invalid or undefined.
        """
        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid.currentText()
        if self.localGridName != "":
            self.layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
            self.diskGridPath = self.layer.dataProvider().dataSourceUri()

            if os.path.exists(self.diskGridPath) or self.localGridName:
                inc, dec, inten = self.getMagParamGeotiff(self.diskGridPath)
                if inc is not None and dec is not None and inten is not None:

                    self.RTE_P_inc = float(inc)
                    self.RTE_P_dec = float(dec)
                    self.RTE_P_int = float(inten)

                    # update widgets
                    self.dlg.lineEdit_5_dec.setText(str(round(self.RTE_P_dec, 1)))
                    self.dlg.lineEdit_6_inc.setText(str(round(self.RTE_P_inc, 1)))
                    self.dlg.lineEdit_6_int.setText(str(int(self.RTE_P_int)))

                else:
                    self.loadGrid()

                    self.base_name = self.localGridName

                    # retrieve parameters
                    self.magn_int = self.dlg.lineEdit_6_int.text()
                    date_text = str(self.dlg.dateEdit.date().toPyDate())

                    date_split = date_text.split("-")
                    self.magn_SurveyDay = int(date_split[2])
                    self.magn_SurveyMonth = int(date_split[1])
                    self.magn_SurveyYear = int(date_split[0])
                    date = datetime(
                        self.magn_SurveyYear, self.magn_SurveyMonth, self.magn_SurveyDay
                    )

                    midx, midy = self.get_grid_centroid(self.layer)

                    if self.layer.isValid():
                        if not self.validCRS(self.layer):
                            return False

                        # convert midpoint to lat/long
                        magn_proj = self.layer.crs().authid().split(":")[1]
                        from pyproj import CRS

                        crs_proj = CRS.from_user_input(int(magn_proj))
                        crs_ll = CRS.from_user_input(4326)
                        proj = Transformer.from_crs(crs_proj, crs_ll, always_xy=True)
                        long, lat = proj.transform(midx, midy)

                        date = self.day_month_to_decimal_year(
                            self.magn_SurveyYear,
                            self.magn_SurveyMonth,
                            self.magn_SurveyDay,
                        )

                        I, D, intensity = self.calcIGRF(date, float(100.0), lat, long)

                        self.RTE_P_inc = I
                        self.RTE_P_dec = D
                        self.RTE_P_int = intensity

                        # update widgets
                        self.dlg.lineEdit_5_dec.setText(str(round(self.RTE_P_dec, 1)))
                        self.dlg.lineEdit_6_inc.setText(str(round(self.RTE_P_inc, 1)))
                        self.dlg.lineEdit_6_int.setText(str(int(self.RTE_P_int)))

    def get_igrf_corners(self):
        """Return (inc_corners, dec_corners, inc_center, dec_center) for the four grid
        corners and the grid centre.

        Corner order matches Cooper_Cowan_rtpvariable:
          [NW, NE, SW, SE]  i.e. [row=nr col=1, row=nr col=nc, row=1 col=1, row=1 col=nc]
          when the grid is displayed with axis-xy (row nr = geographic north).

        Returns (None, None, None, None) on failure.
        """
        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid.currentText()
        if not self.localGridName:
            return None, None, None, None

        self.layer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
        if not self.layer.isValid() or not self.validCRS(self.layer):
            return None, None, None, None

        self.diskGridPath = self.layer.dataProvider().dataSourceUri()
        inc_meta, dec_meta, _ = self.getMagParamGeotiff(self.diskGridPath)

        date_text = str(self.dlg.dateEdit.date().toPyDate())
        date_split = date_text.split("-")
        date = self.day_month_to_decimal_year(
            int(date_split[0]), int(date_split[1]), int(date_split[2])
        )

        if inc_meta is not None:
            # Noddy GeoTIFF: uniform field — centre equals corners
            v_inc, v_dec = float(inc_meta), float(dec_meta)
            return (v_inc,) * 4, (v_dec,) * 4, v_inc, v_dec

        # General case: compute IGRF at each geographic corner and the centre
        magn_proj = self.layer.crs().authid().split(":")[1]
        from pyproj import CRS
        crs_proj = CRS.from_user_input(int(magn_proj))
        crs_ll = CRS.from_user_input(4326)
        proj = Transformer.from_crs(crs_proj, crs_ll, always_xy=True)

        ext = self.layer.extent()
        # Order: [NW, NE, SW, SE] = [xmin/ymax, xmax/ymax, xmin/ymin, xmax/ymin]
        corners_xy = [
            (ext.xMinimum(), ext.yMaximum()),
            (ext.xMaximum(), ext.yMaximum()),
            (ext.xMinimum(), ext.yMinimum()),
            (ext.xMaximum(), ext.yMinimum()),
        ]
        center_xy = ((ext.xMinimum() + ext.xMaximum()) / 2.0,
                     (ext.yMinimum() + ext.yMaximum()) / 2.0)

        inc_corners, dec_corners = [], []
        for x, y in corners_xy:
            lon, lat = proj.transform(x, y)
            inc_val, dec_val, _ = self.calcIGRF(date, 100.0, lat, lon)
            inc_corners.append(inc_val)
            dec_corners.append(dec_val)

        lon_c, lat_c = proj.transform(*center_xy)
        inc_center, dec_center, _ = self.calcIGRF(date, 100.0, lat_c, lon_c)

        if any(v is None for v in inc_corners + dec_corners + [inc_center, dec_center]):
            return None, None, None, None

        print("geographic corners:", corners_xy)
        print("geographic centre:", center_xy)
        print("IGRF corners (inc, dec):", list(zip(inc_corners, dec_corners)))
        print("IGRF centre  (inc, dec):", (inc_center, dec_center))

        return tuple(inc_corners), tuple(dec_corners), inc_center, dec_center

    def getMagParamGeotiff(self, geotiff_path):
        """
        Extract inclination, declination, and intensity from GeoTIFF metadata.

        Parameters:
        geotiff_path (str): Path to the GeoTIFF file

        Returns:
        tuple: (inclination, declination, intensity) as floats, or (None, None, None) if metadata missing

        Raises:
        RuntimeError: If file cannot be opened
        KeyError: If specific metadata keys are missing (unless return_none_if_missing=True)
        """

        # Open the dataset
        dataset = gdal.Open(geotiff_path, gdal.GA_ReadOnly)

        if dataset is None:
            raise RuntimeError(f"Could not open {geotiff_path}")

        # Get metadata dictionary
        metadata = dataset.GetMetadata()

        # Close dataset
        dataset = None

        # Check if required metadata exists
        required_keys = ["inclination", "declination", "intensity"]
        missing_keys = [key for key in required_keys if key not in metadata]

        if missing_keys:
            return None, None, None  # Return None if any key is missing

        # Extract the three parameters and convert to float
        try:
            inclination = float(metadata["inclination"])
            declination = float(metadata["declination"])
            intensity = float(metadata["intensity"])
        except ValueError as e:
            return None, None, None  # Return None if any key is missing

        return inclination, declination, intensity

    def calcIGRF(self, date, alt, lat, lon):
        """
        Calculate the geomagnetic field components and related parameters using the
        International Geomagnetic Reference Field (IGRF) model.

        Parameters:
            date (float): The date for which the geomagnetic field is to be calculated,
                          expressed as a decimal year (e.g., 2023.5 for mid-2023).
            alt (float): The altitude above sea level in kilometers.
            lat (float): The latitude of the location in degrees.
            lon (float): The longitude of the location in degrees.

        Returns:
            tuple: A tuple containing the following:
                - inc (float): Magnetic inclination (angle between the magnetic field
                               vector and the horizontal plane) in degrees.
                - dec (float): Magnetic declination (angle between magnetic north and
                               true north) in degrees.
                - intensity (float): Total intensity of the magnetic field in nanoteslas (nT).

        Notes:
            - The function uses spherical harmonic coefficients from the IGRF model to
              compute the magnetic field components.
            - The coefficients are interpolated to the specified date.
            - The function accounts for secular variation (SV) of the magnetic field
              over time.
            - The calculations include transformations between geocentric and geodetic
              coordinates as needed.
        """

        igrf_gen = "14"
        itype = 1
        d1 = d2 = d3 = None
        colat = 90 - lat
        iut = IGRF(d1, d2, d3)

        # Load in the file of coefficients
        # IGRF_FILE = r"./SHC_files/IGRF" + igrf_gen + ".SHC"
        IGRF_FILE = (
            os.path.dirname(os.path.realpath(__file__))
            + "/calcs/igrf/SHC_files/IGRF"
            + igrf_gen
            + ".SHC"
        )
        from pathlib import Path

        def convert_to_native_path(mixed_path):
            return str(Path(mixed_path))

        IGRF_FILE_norm = convert_to_native_path(IGRF_FILE)
        igrf = iut.load_shcfile(IGRF_FILE_norm, None)

        # Interpolate the geomagnetic coefficients to the desired date(s)
        # -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
        f = interpolate.interp1d(igrf.time, igrf.coeffs, fill_value="extrapolate")
        coeffs = f(date)

        # Compute the main field B_r, B_theta and B_phi value for the location(s)
        Br, Bt, Bp = iut.synth_values(
            coeffs.T, alt, colat, lon, igrf.parameters["nmax"]
        )

        # For the SV, find the 5 year period in which the date lies and compute
        # the SV within that period. IGRF has constant SV between each 5 year period
        # We don't need to subtract 1900 but it makes it clearer:
        epoch = (date - 1900) // 5
        epoch_start = epoch * 5
        # Add 1900 back on plus 1 year to account for SV in nT per year (nT/yr):
        coeffs_sv = f(1900 + epoch_start + 1) - f(1900 + epoch_start)
        Brs, Bts, Bps = iut.synth_values(
            coeffs_sv.T, alt, colat, lon, igrf.parameters["nmax"]
        )

        # Use the main field coefficients from the start of each five epoch
        # to compute the SV for Dec, Inc, Hor and Total Field (F)
        # [Note: these are non-linear components of X, Y and Z so treat separately]
        coeffsm = f(1900 + epoch_start)
        Brm, Btm, Bpm = iut.synth_values(
            coeffsm.T, alt, colat, lon, igrf.parameters["nmax"]
        )

        # Rearrange to X, Y, Z components
        X = -Bt
        Y = Bp
        Z = -Br
        # For the SV
        dX = -Bts
        dY = Bps
        dZ = -Brs
        Xm = -Btm
        Ym = Bpm
        Zm = -Brm
        if itype == 1:
            # alt = input("Enter altitude in km: ").rstrip()
            # alt = iut.check_float(alt)
            alt, colat, sd, cd = iut.gg_to_geo(alt, colat)

        # Rotate back to geodetic coords if needed
        if itype == 1:
            t = X
            X = X * cd + Z * sd
            Z = Z * cd - t * sd
            t = dX
            dX = dX * cd + dZ * sd
            dZ = dZ * cd - t * sd
            t = Xm
            Xm = Xm * cd + Zm * sd
            Zm = Zm * cd - t * sd

        intensity = np.sqrt(X**2 + Y**2 + Z**2)
        # Compute the four non-linear components
        dec, hoz, inc, eff = iut.xyz2dhif(X, Y, Z)
        return inc, dec, intensity

    def extract_raster_to_numpy(self, raster_layer):
        """
        Extract the raster data from a QgsRasterLayer as a NumPy array.

        Parameters:
            raster_layer (QgsRasterLayer): The QGIS raster layer.

        Returns:
            numpy.ndarray: The raster data as a NumPy array.
        """
        # Get the raster data provider
        """
        # Access the raster data provider
        provider = raster_layer.dataProvider()

        # Get raster dimensions
        cols = provider.xSize()  # Number of columns
        rows = provider.ySize()  # Number of rows

        # Read raster data as a block
        band = 1  # Specify the band number (1-based index)
        raster_block = provider.block(band, provider.extent(), cols, rows)

        # Copy the block data into a NumPy array
        extent = raster_layer.extent()
        rows, cols = raster_layer.height(), raster_layer.width()
        raster_block = provider.block(1, extent, cols, rows)  # !!!!!
        self.raster_array = np.zeros((rows, cols))
        for i in range(rows):
            for j in range(cols):
                self.raster_array[i, j] = raster_block.value(i, j)
        """
        # Read raster data via GDAL (consistent row ordering across platforms)
        path = raster_layer.dataProvider().dataSourceUri().split("|")[0]
        ds = gdal.Open(path)
        band = ds.GetRasterBand(1)
        no_data_value = band.GetNoDataValue()  # Band 1
        band_data = ds.GetRasterBand(1).ReadAsArray()  # Always row0=north
        if no_data_value is not None:
            band_data = band_data.astype(float)
            band_data[band_data == no_data_value] = np.nan
        self.raster_array = band_data
        gt = ds.GetGeoTransform()
        if gt[5] > 0:  # positive pixel height means south-up, flip to north-up
            self.raster_array = np.flipud(self.raster_array)
        ds = None

        return self.raster_array

    def display_rad_power_spectrum(self):
        """
        Displays the radial power spectrum of a selected raster layer in QGIS.
        This method performs the following steps:
        1. Checks for the required Python packages (`pywt` and `matplotlib`) and prompts the user to install them if missing.
        2. Retrieves the selected raster layer from the QGIS interface.
        3. Extracts the raster data as a NumPy array.
        4. Computes the spatial resolution and extent of the raster layer.
        5. Initializes a PowerSpectrumDock instance to visualize the raster grid and its power spectrum.
        Returns:
            bool: False if required packages are missing or if the selected raster layer is invalid.
        Raises:
            None: This method handles exceptions related to missing packages and invalid layers internally.
        Notes:
            - The method assumes the existence of a method `extract_raster_to_numpy` to convert raster data to a NumPy array.
            - The `PowerSpectrumDock` class is used for plotting the grid and power spectrum.
            - This method interacts with the QGIS Python API and requires a valid QGIS environment.
        Dependencies:
            - pywt: For wavelet transformations.
            - matplotlib: For plotting.
        Example:
            To use this method, ensure that a raster layer is selected in the QGIS interface and the required packages are installed.
        """

        try:
            import matplotlib
        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for some functions, but not installed: matplotlib\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install matplotlib",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False

        try:
            import matplotlib.colors as mcolors

        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for some functions, but not installed: matplotlib\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install matplotlib",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False

        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid.currentText()
        if self.localGridName != "":
            self.pslayer = QgsProject.instance().mapLayersByName(self.localGridName)[0]
            if self.pslayer.isValid():
                grid = self.extract_raster_to_numpy(
                    self.pslayer
                )  # Your method to get NumPy array from raster

                dx, dy = (
                    self.pslayer.rasterUnitsPerPixelX(),
                    self.pslayer.rasterUnitsPerPixelY(),
                )
                # Get extent
                extent = self.pslayer.extent()
                minx = extent.xMinimum()
                maxx = extent.xMaximum()
                miny = extent.yMinimum()
                maxy = extent.yMaximum()

                # Get number of columns (nx) and rows (ny)
                provider = self.pslayer.dataProvider()
                nx = provider.xSize()  # Number of columns
                ny = provider.ySize()  # Number of rows

                x = np.linspace(minx, maxx, provider.xSize())
                y = np.linspace(miny, maxy, provider.ySize())
                # Initialize the PowerSpectrumDock and display the plot
                power_spectrum_dock = PowerSpectrumDock(
                    grid, self.localGridName, dx, dy, x, y
                )
                power_spectrum_dock.plot_grid_and_power_spectrum()

    def update_paths(self):
        """
        Updates the paths and UI elements related to the selected grid layer.
        This method performs the following actions:
        - Retrieves the name of the currently selected grid layer from the UI.
        - Updates multiple combo boxes in the UI to reflect the selected grid layer.
        - Clears the grid path line edit field and resets the disk grid path.
        - Sets the base name of the grid layer to the selected grid name.
        - Checks if the selected grid layer exists and is valid.
        - Determines the coordinate reference system (CRS) of the selected layer.
        - Updates the units label in the UI based on whether the CRS is geographic
          (degrees) or projected (meters).
        Note:
            This method assumes that the selected grid layer exists in the current
            QGIS project and that the UI elements are properly initialized.
        Raises:
            IndexError: If no layer with the selected grid name is found in the
                        current QGIS project.
        """

        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid.currentText()
        self.dlg.mMapLayerComboBox_selectGrid_Conv.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_worms.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_Conv_2.setCurrentText(self.localGridName)
        self.dlg.lineEdit_2_loadGridPath.setText("")
        self.diskGridPath = ""
        self.base_name = self.localGridName

        if len(self.base_name) > 0:
            selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[
                0
            ]
            if selected_layer.isValid():
                crs = selected_layer.crs()
                if crs.isGeographic():
                    self.dlg.label_41_units.setText("Units: deg")
                else:
                    self.dlg.label_41_units.setText("Units: m")

    def update_paths_utils(self):
        """
        Updates the paths and UI elements related to the selected grid layer.
        This method performs the following actions:
        - Retrieves the name of the currently selected grid layer from the UI.
        - Updates multiple combo boxes in the UI to reflect the selected grid layer.
        - Clears the grid path line edit field and resets the disk grid path.
        - Sets the base name of the grid layer to the selected grid name.
        - Checks if the selected grid layer exists and is valid.
        - Determines the coordinate reference system (CRS) of the selected layer.
        - Updates the units label in the UI based on whether the CRS is geographic
          (degrees) or projected (meters).
        Note:
            This method assumes that the selected grid layer exists in the current
            QGIS project and that the UI elements are properly initialized.
        Raises:
            IndexError: If no layer with the selected grid name is found in the
                        current QGIS project.
        """
        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid_Conv_2.currentText()
        self.dlg.mMapLayerComboBox_selectGrid_Conv.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_worms.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid.setCurrentText(self.localGridName)

        self.dlg.lineEdit_2_loadGridPath.setText("")
        self.diskGridPath = ""
        self.base_name = self.localGridName

        if len(self.base_name) > 0:
            selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[
                0
            ]
            if selected_layer.isValid():
                crs = selected_layer.crs()
                if crs.isGeographic():
                    self.dlg.label_41_units.setText("Units: deg")
                else:
                    self.dlg.label_41_units.setText("Units: m")

    def update_paths_conv(self):
        """
        Updates the paths and UI elements related to the selected grid layer.
        This method performs the following actions:
        - Retrieves the name of the currently selected grid layer from the UI.
        - Updates multiple combo boxes in the UI to reflect the selected grid layer.
        - Clears the grid path line edit field and resets the disk grid path.
        - Sets the base name of the grid layer to the selected grid name.
        - Checks if the selected grid layer exists and is valid.
        - Determines the coordinate reference system (CRS) of the selected layer.
        - Updates the units label in the UI based on whether the CRS is geographic
          (degrees) or projected (meters).
        Note:
            This method assumes that the selected grid layer exists in the current
            QGIS project and that the UI elements are properly initialized.
        Raises:
            IndexError: If no layer with the selected grid name is found in the
                        current QGIS project.
        """
        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid_Conv.currentText()
        self.dlg.mMapLayerComboBox_selectGrid.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_worms.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_Conv_2.setCurrentText(self.localGridName)

        self.dlg.lineEdit_2_loadGridPath.setText("")
        self.diskGridPath = ""
        self.base_name = self.localGridName

        if len(self.base_name) > 0:
            selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[
                0
            ]
            if selected_layer.isValid():
                crs = selected_layer.crs()
                if crs.isGeographic():
                    self.dlg.label_41_units.setText("Units: deg")
                else:
                    self.dlg.label_41_units.setText("Units: m")

    def update_paths_worms(self):
        """
        Updates the paths and UI elements related to the selected grid layer.
        This method performs the following actions:
        - Retrieves the name of the currently selected grid layer from the UI.
        - Updates multiple combo boxes in the UI to reflect the selected grid layer.
        - Clears the grid path line edit field and resets the disk grid path.
        - Sets the base name of the grid layer to the selected grid name.
        - Checks if the selected grid layer exists and is valid.
        - Determines the coordinate reference system (CRS) of the selected layer.
        - Updates the units label in the UI based on whether the CRS is geographic
          (degrees) or projected (meters).
        Note:
            This method assumes that the selected grid layer exists in the current
            QGIS project and that the UI elements are properly initialized.
        Raises:
            IndexError: If no layer with the selected grid name is found in the
                        current QGIS project.
        """
        self.localGridName = self.dlg.mMapLayerComboBox_selectGrid_worms.currentText()

        self.dlg.mMapLayerComboBox_selectGrid.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_Conv.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_Conv_2.setCurrentText(self.localGridName)
        self.dlg.mMapLayerComboBox_selectGrid_worms.setCurrentText(self.localGridName)

        self.dlg.lineEdit_2_loadGridPath.setText("")
        self.diskGridPath = ""
        self.base_name = self.localGridName

        if len(self.base_name) > 0:
            selected_layer = QgsProject.instance().mapLayersByName(self.localGridName)[
                0
            ]
            if selected_layer.isValid():
                crs = selected_layer.crs()
                if crs.isGeographic():
                    self.dlg.label_41_units.setText("Units: deg")
                else:
                    self.dlg.label_41_units.setText("Units: m")

    # --------------------------------------------------------------------------
    def show_version(self):
        """Show the version of the plugin"""
        metadata_path = os.path.dirname(os.path.realpath(__file__)) + "/metadata.txt"

        with open(metadata_path) as plugin_version_file:
            metadata = plugin_version_file.readlines()
            for line in metadata:
                parts = line.split("=")
                if len(parts) == 2 and parts[0] == "version":
                    plugin_version = parts[1]

            return plugin_version

    def run(self):
        """Run method that loads and starts the plugin"""

        if not self.pluginIsActive:
            self.pluginIsActive = True
            self.initParams()
            # print "** STARTING SGTool"

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

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

            # show the dockwidget
            # TODO: fix to allow choice of dock location
            # QT5 self.iface.addDockWidget(Qt.RightDockWidgetArea, self.dlg)
            self.iface.addDockWidget(RightDockWidgetArea, self.dlg)
            # 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.dlg:
                        self.iface.mainWindow().tabifyDockWidget(dock, self.dlg)
                        # Optionally, bring your plugin tab to the front
                        self.dlg.raise_()
                        break
            # Raise the docked widget above others
            self.dlg.show()
            self.define_tips()

            # Access the QgsMapLayerComboBox by its objectName
            self.dlg.mMapLayerComboBox_selectGrid.setFilters(
                QgsMapLayerProxyModel.RasterLayer
            )
            self.dlg.mMapLayerComboBox_selectGrid_Conv.setFilters(
                QgsMapLayerProxyModel.RasterLayer
            )
            self.dlg.mMapLayerComboBox_selectGrid_worms.setFilters(
                QgsMapLayerProxyModel.RasterLayer
            )
            self.dlg.mMapLayerComboBox_selectGrid_Conv_2.setFilters(
                QgsMapLayerProxyModel.RasterLayer
            )

            self.dlg.mMapLayerComboBox_selectVectors.setFilters(
                QgsMapLayerProxyModel.PointLayer | QgsMapLayerProxyModel.VectorLayer
            )

            self.dlg.mMapLayerComboBox_selectGrid_3.setFilters(
                QgsMapLayerProxyModel.PointLayer
            )

            self.dlg.version_label.setText(self.show_version())

            self.deriv_dir_list = []
            self.deriv_dir_list.append("z")
            self.deriv_dir_list.append("x")
            self.deriv_dir_list.append("y")
            self.dlg.comboBox_derivDirection.addItems(self.deriv_dir_list)

            self.contin_dir_list = []
            self.contin_dir_list.append("up")
            self.contin_dir_list.append("down")
            self.dlg.comboBox_2_continuationDirection.addItems(self.contin_dir_list)

            self.dlg.comboBox_3_rte_p_list.clear()
            self.ret_p_list = ["RTP", "Diff. RTP", "RTE"]
            self.dlg.comboBox_3_rte_p_list.addItems(self.ret_p_list)

            self.freq_cut_type_list = []
            self.freq_cut_type_list.append("Low")
            self.freq_cut_type_list.append("High")
            self.dlg.comboBox_2_FreqCutType.addItems(self.freq_cut_type_list)

            self.dlg.pushButton_4_calcIGRF.clicked.connect(self.update_mag_field)
            self.dlg.pushButton_2_selectGrid.clicked.connect(self.select_grid_file)
            self.dlg.pushButton_2_selectGrid_RGB.clicked.connect(
                self.select_RGBgrid_file
            )
            self.dlg.pushButton_3_applyProcessing.clicked.connect(
                self.processGeophysics_fft
            )
            self.update_wavelet_choices()

            self.dlg.pushButton_3_applyProcessing_Conv_3.clicked.connect(
                self.processGeophysics_fft
            )
            self.dlg.pushButton_3_applyProcessing_Conv.clicked.connect(
                self.processGeophysics_conv
            )
            self.dlg.pushButton_3_applyProcessing_Conv_2.clicked.connect(
                self.processRGB
            )

            self.dlg.pushButton_selectPoints.clicked.connect(self.select_point_file)

            self.dlg.mMapLayerComboBox_selectGrid.layerChanged.connect(
                self.update_paths
            )
            self.dlg.mMapLayerComboBox_selectGrid_Conv.layerChanged.connect(
                self.update_paths_conv
            )
            self.dlg.mMapLayerComboBox_selectGrid_Conv_2.layerChanged.connect(
                self.update_paths_utils
            )

            self.dlg.mMapLayerComboBox_selectVectors.layerChanged.connect(
                self.update_wavelet_choices
            )
            self.dlg.mMapLayerComboBox_selectGrid_worms.layerChanged.connect(
                self.update_paths_worms
            )

            self.dlg.mQgsProjectionSelectionWidget.setCrs(
                QgsCoordinateReferenceSystem("EPSG:4326")
            )
            self.dlg.pushButton_rad_power_spectrum.clicked.connect(
                self.display_rad_power_spectrum
            )
            self.localGridName = self.dlg.mMapLayerComboBox_selectGrid.currentText()

            if self.localGridName:
                selected_layer = QgsProject.instance().mapLayersByName(
                    self.localGridName
                )
                if not selected_layer:
                    return

                selected_layer = selected_layer[0]
                crs = selected_layer.crs()
                if crs.isGeographic():
                    self.dlg.label_41_units.setText("Units: deg")
                else:
                    self.dlg.label_41_units.setText("Units: m")

            self.dlg.lineEdit_Mean_size.setValidator(OddPositiveIntegerValidator())
            self.dlg.lineEdit_Median_size.setValidator(OddPositiveIntegerValidator())
            self.dlg.lineEdit_SS_Window.setValidator(OddPositiveIntegerValidator())
            self.dlg.lineEdit_ED_Window.setValidator(OddPositiveIntegerValidator())

            self.directional_list = []
            self.directional_list.append("N")
            self.directional_list.append("NE")
            self.directional_list.append("E")
            self.directional_list.append("SE")
            self.directional_list.append("S")
            self.directional_list.append("SW")
            self.directional_list.append("W")
            self.directional_list.append("NW")
            self.dlg.comboBox_Dir_dir.addItems(self.directional_list)

            self.dlg.pushButton_load_point_data.clicked.connect(
                self.import_point_line_data
            )
            self.dlg.mMapLayerComboBox_selectGrid_3.layerChanged.connect(
                self.updateLayertoGrid
            )

            self.dlg.doubleSpinBox_cellsize.valueChanged.connect(
                self.updateLayertoGrid2
            )
            self.dlg.comboBox_select_grid_data_field.currentTextChanged.connect(
                self.updateLayertoGrid2
            )

            if self.dlg.mMapLayerComboBox_selectGrid_3.currentText() != "":
                self.updateLayertoGrid()

            self.cell_size = self.dlg.doubleSpinBox_cellsize.value()

            self.gridDirectory = None

            # Connection to the CSS Colours site  :
            self.dlg.pushButton_CSSS_Colours.clicked.connect(
                lambda: QDesktopServices.openUrl(
                    QUrl(
                        "https://matplotlib.org/stable/gallery/color/named_colors.html#css-colors"
                    )
                )
            )

            self.dlg.pushButton_pick_rgb.clicked.connect(self._on_pick_rgb_clicked)

            self.dlg.pushButton_idw_2.clicked.connect(self.procIDWGridding)
            self.dlg.pushButton_bspline_3.clicked.connect(self.procmultibsplineGridding)

            self.dlg.pushButton_worms.clicked.connect(self.procBSDworms)

            self.dlg.pushButton_normalise.clicked.connect(self.procNormalise)
            self.dlg.pushButton_select_normalise_in.clicked.connect(
                self.set_normalise_in
            )
            self.dlg.pushButton_select_normalise_out.clicked.connect(
                self.set_normalise_out
            )

            self.dlg.pushButton_wtmm.clicked.connect(self.procWTMM)
            # autocheck the associated checkbox if a parameter is modified

            self.dlg.lineEdit_3_azimuth.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_3_DirClean)
            )
            self.dlg.lineEdit_3_DC_wavelength.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_3_DirClean)
            )
            self.dlg.lineEdit_3_DC_scale.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_3_DirClean)
            )
            self.dlg.lineEdit_6_int.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_4_RTE_P)
            )
            self.dlg.lineEdit_6_inc.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_4_RTE_P)
            )
            self.dlg.lineEdit_5_dec.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_4_RTE_P)
            )
            self.dlg.dateEdit.dateChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_4_RTE_P)
            )
            self.dlg.comboBox_3_rte_p_list.currentTextChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_4_RTE_P)
            )
            self.dlg.comboBox_2_continuationDirection.currentTextChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_9_continuation)
            )
            self.dlg.lineEdit_10_continuationHeight.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_9_continuation)
            )
            self.dlg.radioButton_RR_1st.clicked.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_5_regional)
            )
            self.dlg.radioButton_RR_2nd.clicked.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_5_regional)
            )
            self.dlg.lineEdit_12_bandPassLow.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_10_bandPass)
            )
            self.dlg.lineEdit_11_bandPassHigh.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_10_bandPass)
            )
            self.dlg.lineEdit_3_BP_width.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_10_bandPass)
            )
            self.dlg.lineEdit_3_HLP_width.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_10_freqCut)
            )
            self.dlg.lineEdit_12_FreqPass.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_10_freqCut)
            )
            self.dlg.comboBox_2_FreqCutType.currentTextChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_10_freqCut)
            )
            self.dlg.lineEdit_13_agc_window.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_11_1vd_agc)
            )
            self.dlg.lineEdit_9_derivePower.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_6_derivative)
            )
            self.dlg.comboBox_derivDirection.currentTextChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_6_derivative)
            )
            self.dlg.lineEdit_Mean_size.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_Mean)
            )
            self.dlg.lineEdit_Median_size.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_Median)
            )
            self.dlg.lineEdit_Gaussian_Sigma.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_Gaussian)
            )
            self.dlg.comboBox_Dir_dir.currentTextChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_Directional)
            )
            self.dlg.checkBox_relief.toggled.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_SunShading)
            )
            self.dlg.lineEdit_SunSh_Az.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_SunShading)
            )
            self.dlg.lineEdit_SunSh_Zn.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_SunShading)
            )
            self.dlg.radioButton_NaN_Above.toggled.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_NaN)
            )
            self.dlg.radioButton_NaN_Below.toggled.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_NaN)
            )
            self.dlg.radioButton_NaN_Both.toggled.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_NaN)
            )
            self.dlg.doubleSpinBox_NaN_Above.valueChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_NaN)
            )
            self.dlg.doubleSpinBox_NaN_Below.valueChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_NaN)
            )
            for _mrvbf_edit in (
                self.dlg.lineEdit_MRVBF_tslope,
                self.dlg.lineEdit_MRVBF_tpctlv,
                self.dlg.lineEdit_MRVBF_tpctlr,
                self.dlg.lineEdit_MRVBF_pslope,
                self.dlg.lineEdit_MRVBF_ppctl,
                self.dlg.lineEdit_MRVBF_maxres,
            ):
                _mrvbf_edit.textChanged.connect(
                    lambda: self.update_checkbox(self.dlg.checkBox_MRVBF)
                )
            self.dlg.mQgsSpinBox_PCA.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_PCA)
            )
            self.dlg.mQgsSpinBox_ICA.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_ICA)
            )
            self.dlg.checkBox_ED_Stats.toggled.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_ED)
            )
            self.dlg.doubleSpinBox_ED_Threshold.valueChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_ED)
            )
            self.dlg.lineEdit_ED_Window.textChanged.connect(
                lambda: self.update_checkbox(self.dlg.checkBox_ED)
            )

    def update_checkbox(self, checkbox):
        """Generic method to set a checkbox to checked state."""
        checkbox.setChecked(True)

    def get_layer_fields(self, layer):
        """
        Get a list of field names from a QgsVectorLayer.

        Parameters:
            layer (QgsVectorLayer): The vector layer to retrieve fields from.

        Returns:
            list: A list of field names.
        """
        if not layer.isValid():
            raise ValueError("Invalid layer provided.")

        fields = layer.fields()
        field_names = [field.name() for field in fields]
        return field_names

    def updateLayertoGrid(self):
        """
        Updates the grid layer selection and calculates grid dimensions based on the selected layer's extent
        and the specified cell size.
        This method performs the following steps:
        1. Checks if there are any layers available in the grid selection combo box.
        2. Retrieves the currently selected layer from the combo box.
        3. Validates the selected layer and retrieves its field names to populate the grid data field combo box.
        4. If the selected layer contains features, calculates the grid dimensions (nx_label and ny_label)
           based on the layer's extent and the user-defined cell size.
        5. Updates the corresponding labels in the dialog with the calculated grid dimensions.
        Attributes:
            dlg.mMapLayerComboBox_selectGrid_3 (QComboBox): Combo box for selecting the grid layer.
            dlg.comboBox_select_grid_data_field (QComboBox): Combo box for selecting grid data fields.
            dlg.doubleSpinBox_cellsize (QDoubleSpinBox): Spin box for specifying the cell size.
            dlg.nx_label (QLabel): Label to display the calculated number of grid cells along the x-axis.
            dlg.ny_label (QLabel): Label to display the calculated number of grid cells along the y-axis.
        Raises:
            IndexError: If no layers are found with the selected name in the QgsProject instance.
            AttributeError: If the selected layer is invalid or does not have the required attributes.
        """

        if self.dlg.mMapLayerComboBox_selectGrid_3.count() > 0:

            self.selectedPoints = self.dlg.mMapLayerComboBox_selectGrid_3.currentText()

            layers = QgsProject.instance().mapLayersByName(self.selectedPoints)

            if not layers:
                return

            selected_layer = layers[0]

            if selected_layer.isValid():
                field_names = self.get_layer_fields(selected_layer)
                self.dlg.comboBox_select_grid_data_field.clear()
                self.dlg.comboBox_select_grid_data_field.addItems(field_names)

                if selected_layer.featureCount() > 0:
                    extent = selected_layer.extent()

                    self.cell_size = self.dlg.doubleSpinBox_cellsize.value()

                    _, _, self.nx_label, self.ny_label = self._mba_grid_geometry(
                        extent.xMinimum(), extent.xMaximum(),
                        extent.yMinimum(), extent.yMaximum(), self.cell_size,
                    )
                    self.dlg.nx_label.setText(str(self.nx_label))
                    self.dlg.ny_label.setText(str(self.ny_label))

    def updateLayertoGrid2(self):
        """
        Updates the grid dimensions (nx_label and ny_label) based on the selected layer
        and cell size from the user interface.

        This method retrieves the currently selected layer from the
        `mMapLayerComboBox_selectGrid_3` combo box, calculates the grid dimensions
        (number of cells along the x and y axes) based on the layer's extent and
        the specified cell size, and updates the corresponding labels in the dialog.

        Preconditions:
            - The combo box `mMapLayerComboBox_selectGrid_3` must contain at least one layer.
            - The selected layer must be valid and contain features.

        Steps:
            1. Retrieve the selected layer from the combo box.
            2. Check if the layer is valid.
            3. Calculate the grid dimensions (nx_label and ny_label) using the layer's extent
               and the cell size specified in `doubleSpinBox_cellsize`.
            4. Update the `nx_label` and `ny_label` text fields in the dialog.

        Attributes:
            dlg.mMapLayerComboBox_selectGrid_3 (QComboBox): Combo box for selecting a grid layer.
            dlg.doubleSpinBox_cellsize (QDoubleSpinBox): Spin box for specifying the cell size.
            dlg.nx_label (QLabel): Label to display the number of cells along the x-axis.
            dlg.ny_label (QLabel): Label to display the number of cells along the y-axis.

        Raises:
            IndexError: If no layers are found with the selected name.
            AttributeError: If the selected layer is invalid or does not have the required attributes.

        """

        if self.dlg.mMapLayerComboBox_selectGrid_3.count() > 0:

            self.selectedPoints = self.dlg.mMapLayerComboBox_selectGrid_3.currentText()
            selected_layer = QgsProject.instance().mapLayersByName(self.selectedPoints)
            if not selected_layer:
                return

            selected_layer = selected_layer[0]

            if selected_layer.isValid():

                extent = selected_layer.extent()
                self.cell_size = self.dlg.doubleSpinBox_cellsize.value()

                if selected_layer.featureCount() > 0:
                    _, _, self.nx_label, self.ny_label = self._mba_grid_geometry(
                        extent.xMinimum(), extent.xMaximum(),
                        extent.yMinimum(), extent.yMaximum(), self.cell_size,
                    )
                    self.dlg.nx_label.setText(str(self.nx_label))
                    self.dlg.ny_label.setText(str(self.ny_label))

    def import_point_line_data(self, header_list=None, data=None):
        # import point or line data as vector file to memory
        self.diskPointsPath = self.dlg.lineEdit_loadPointsPath.text()
        if os.path.exists(self.diskPointsPath) and self.diskPointsPath != "":

            dir_name, base_name = os.path.split(self.diskPointsPath)
            file_name, file_ext = os.path.splitext(base_name)

            proj = self.dlg.mQgsProjectionSelectionWidget.crs().authid()

            if self.pointType == "line":
                crs = proj.split(":")[1]
                load_ties = self.dlg.checkBox_load_tie_lines.isChecked()
                self.import_XYZ(
                    self.diskPointsPath, crs, file_name, load_ties=load_ties
                )
            else:
                x_field = self.dlg.comboBox_grid_x.currentText()
                y_field = self.dlg.comboBox_grid_y.currentText()
                if file_ext.upper() == ".CSV" or file_ext.upper() == ".TXT":
                    self.import_CSV(
                        self.diskPointsPath,
                        x_field,
                        y_field,
                        layer_name=file_name,
                        crs=proj,
                    )
                elif file_ext.upper() == ".DAT":
                    data = self.parser.parse_dat_file(
                        self.diskPointsPath, self.header_list, self.field_defs
                    )

                    if not self.points_epsg:
                        self.points_epsg = proj.split(":")[1]

                    if data:
                        # Check if the layer exists, delete if it does
                        previous_layer = QgsProject.instance().mapLayersByName(
                            file_name
                        )
                        print("previous_layer", previous_layer)
                        if len(previous_layer) > 0:
                            previous_layer = previous_layer[0]
                            if previous_layer.isValid():
                                # Store the layer ID
                                layer_id = previous_layer.id()

                                # Remove the layer from registry but don't delete it yet
                                QgsProject.instance().removeMapLayer(layer_id)
                                print("previous layer id removed", layer_id)

                                # Force garbage collection to release file locks
                                import gc

                                previous_layer = None  # Remove reference to the layer
                                gc.collect()  # Force garbage collection
                                import time

                                time.sleep(0.5)  # Wait half a second

                        output_file = os.path.join(dir_name, file_name + ".shp")

                        output_path = self.parser.export_to_shapefile(
                            data,
                            self.field_defs,
                            output_file,
                            file_name,
                            x_field,
                            y_field,
                            self.points_epsg,
                        )

                        layer = QgsVectorLayer(output_path, file_name)
                        # Add the layer to the current QGIS project
                        QgsProject.instance().addMapLayer(layer)

                        canvas = (
                            self.iface.mapCanvas() if hasattr(self, "iface") else None
                        )
                        extent = layer.extent()
                        buffer_amount = max(extent.width(), extent.height()) * 0.05
                        buffered_extent = extent.buffered(buffer_amount)
                        canvas.setExtent(buffered_extent)
                        canvas.refresh()

                    else:
                        print("Error: No data found in the file.")

                    """self.create_points_layer_from_data(
                        dir_name,
                        self.pts_columns,
                        proj.split(":")[1],
                        self.pts_data,
                        x_field,
                        y_field,
                        layer_name=file_name,
                    )"""

    def import_CSV(
        self, file_path, x_field, y_field, layer_name="points", crs="EPSG:4326"
    ):
        """
        Loads a CSV file as a vector layer with all attributes in QGIS.

        Parameters:
            file_path (str): Path to the CSV file.
            x_field (str): Name of the column containing X coordinates.
            y_field (str): Name of the column containing Y coordinates.
            layer_name (str): Name for the layer in the QGIS project.
            crs (str): Coordinate reference system for the layer (default is 'EPSG:4326').

        Returns:
            QgsVectorLayer: The loaded vector layer.
        """
        # Define the URI for the CSV file, specifying coordinate fields and CRS

        dir_name, basename = os.path.split(file_path)
        output_path = os.path.join(dir_name, f"{basename}.shp")

        uri = (
            f"file:///{file_path}?type=csv&xField={x_field}&yField={y_field}"
            f"&crs={crs}&detectTypes=yes&delimiter=,&quote="
        )

        # Load the layer as a delimited text vector layer
        layer = QgsVectorLayer(uri, layer_name, "delimitedtext")
        if not layer.isValid():
            raise ValueError(f"Failed to load layer: {file_path}")
        self.save_layer_as_shapefile(layer, output_path)

        layer2 = QgsVectorLayer(output_path, layer_name)
        # Add the layer to the current QGIS project
        QgsProject.instance().addMapLayer(layer2)
        return layer2

    def save_layer_as_shapefile(self, layer, output_path):
        """
        Save a vector layer as a shapefile using QgsVectorFileWriter3.

        Parameters:
        layer -- The QgsVectorLayer to save
        output_path -- The full path where the shapefile will be saved (including .shp extension)

        Returns:
        bool -- True if successful, False otherwise
        """
        # Make sure the layer is valid
        if not layer.isValid():
            print("Layer is not valid")
            return False

        # Create the save options
        options = QgsVectorFileWriter.SaveVectorOptions()
        options.driverName = "ESRI Shapefile"
        options.fileEncoding = "UTF-8"

        # Save the layer using QgsVectorFileWriter3
        write_result, error_message, new_filename, new_layer_name = (
            QgsVectorFileWriter.writeAsVectorFormatV3(
                layer, output_path, QgsCoordinateTransformContext(), options
            )
        )
        # Check if the saving was successful
        if write_result == QgsVectorFileWriter.NoError:
            print(f"Layer successfully saved to {output_path}")
            return True
        else:
            print(f"Error saving layer: {error_message}")
            return False

    def create_points_layer_from_data(
        self,
        dir_name,
        header_list,
        epsg,
        data,
        x_field_name=None,
        y_field_name=None,
        layer_name="Points Layer",
    ):
        """
        Create a shapefile from the data dictionary and load it as a layer.
        """

        # Determine coordinate field names if not provided
        if x_field_name is None:
            # Try to find X field in header list
            for field in header_list:
                if "X_" in field or field.endswith("_X") or "EAST" in field.upper():
                    x_field_name = field
                    print(f"Using {x_field_name} as X coordinate field")
                    break

        if y_field_name is None:
            # Try to find Y field in header list
            for field in header_list:
                if "Y_" in field or field.endswith("_Y") or "NORTH" in field.upper():
                    y_field_name = field
                    print(f"Using {y_field_name} as Y coordinate field")
                    break

        # Verify we have coordinate fields
        if x_field_name is None or y_field_name is None:
            self.iface.messageBar().pushMessage(
                "Could not determine coordinate fields", level=Qgis.Warning, duration=3
            )
            return None

        # Create a CRS object
        crs = QgsCoordinateReferenceSystem("EPSG:4326")  # Default to WGS84
        if epsg is not None:
            crs = QgsCoordinateReferenceSystem(f"EPSG:{epsg}")
            print(f"Using CRS: {crs.description()}")

        # Create fields
        fields = QgsFields()
        for field_name in header_list:
            # Determine field type from data
            field_type = QVariant.String  # Default to string

            for row in data:
                if field_name in row and row[field_name] is not None:
                    value = row[field_name]
                    if isinstance(value, int):
                        field_type = QVariant.Int
                        break
                    elif isinstance(value, float):
                        field_type = QVariant.Double
                        break

            fields.append(QgsField(field_name, field_type))

        # Create a temporary shapefile
        # temp_dir = tempfile.gettempdir()
        output_file = os.path.join(dir_name, f"{layer_name}.shp")

        print(f"Creating shapefile: {output_file}")

        # Create the shapefile writer
        writer_options = QgsVectorFileWriter.SaveVectorOptions()
        writer_options.driverName = "ESRI Shapefile"
        writer_options.fileEncoding = "UTF-8"

        transform_context = QgsProject.instance().transformContext()

        writer = QgsVectorFileWriter.create(
            output_file,
            fields,
            QgsWkbTypes.Point,
            crs,
            transform_context,
            writer_options,
        )

        # Add features
        feature_count = 0
        error_count = 0

        # Verify that coordinate fields exist in the data
        if data and (x_field_name not in data[0] or y_field_name not in data[0]):
            print(
                f"Warning: coordinate fields not found in data. Available fields: {list(data[0].keys())}"
            )
            return None

        for i, row in enumerate(data):
            try:
                # Get coordinates
                x_coord = float(row[x_field_name])
                y_coord = float(row[y_field_name])

                # Create feature
                feat = QgsFeature(fields)

                # Set geometry
                feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(x_coord, y_coord)))

                # Set attributes
                for j, field_name in enumerate(header_list):
                    feat.setAttribute(j, row.get(field_name, None))

                # Add to shapefile
                writer.addFeature(feat)
                feature_count += 1

            except Exception as e:
                error_count += 1
                if error_count < 10:  # Limit the number of error messages
                    print(f"Error with row {i}: {e}")

        # Clean up writer
        del writer

        print(f"Features written to shapefile: {feature_count}")
        print(f"Errors encountered: {error_count}")

        # Load the shapefile
        layer = QgsVectorLayer(output_file, layer_name, "ogr")

        if not layer.isValid():
            print(f"Layer is not valid! Path: {output_file}")
            return None

        print(f"Shapefile loaded with {layer.featureCount()} features")

        # Add to project
        QgsProject.instance().addMapLayer(layer)

        # Try to zoom to layer
        try:
            self.iface.setActiveLayer(layer)
            self.iface.zoomToActiveLayer()
        except Exception as e:
            print(f"Could not zoom to layer: {str(e)}")

        return layer

    def import_XYZ(self, XYZ_file, crs, layer_name="line", load_ties=True):
        """
        Imports an XYZ file and creates both line and point layers in QGIS.
        This function reads an XYZ file containing coordinate data and optional line or tie markers.
        It processes the data to create a line layer and a point layer in memory, which are added
        to the QGIS project. The line layer represents connected points for each line ID, while the
        point layer represents individual points.
        Args:
            XYZ_file (str): The file path to the XYZ file to be imported.
            crs (str): The coordinate reference system (CRS) in EPSG code format (e.g., "4326").
            layer_name (str, optional): The base name for the created layers. Defaults to "line".
            load_ties (bool, optional): Whether to load tie points (lines marked with "TIE:").
                                        Defaults to True.
        Raises:
            ValueError: If the XYZ file contains invalid data that cannot be parsed.
        Notes:
            - The XYZ file should contain lines with coordinates (x, y, [optional additional data])
              and optional markers "LINE:" or "TIE:" to indicate line IDs.
            - The created layers are added to the current QGIS project.
            - The point layer is initially set to be invisible in the layer tree.
        Example:
            import_XYZ("path/to/xyz_file.xyz", "4326", layer_name="my_layer", load_ties=False)
        """

        # Initialize variables
        data_list = []
        current_line_number = None
        dir_name, basename = os.path.split(XYZ_file)
        output_path = os.path.join(dir_name, f"{basename}.shp")

        # Read the file line-by-line
        with open(XYZ_file, "r") as file:
            for line in file:
                line = line.strip()
                if line.startswith("LINE:"):  # Check for 'LINE:' markers
                    current_line_number = int(re.search(r"\d+", line).group())
                elif line.startswith("TIE:"):  # Check for 'TIE:' markers
                    if load_ties:
                        current_line_number = int(re.search(r"\d+", line).group())
                    else:
                        current_line_number = None
                elif current_line_number is not None:
                    try:
                        parts = list(map(float, line.split()))
                        if len(parts) >= 2:  # Ensure at least x and y are present
                            data_list.append(parts + [current_line_number])
                    except ValueError:
                        pass
                else:
                    if load_ties:
                        print("Invalid line:", line)
        # Process and create the line layer
        line_layer = QgsVectorLayer("LineString?crs=EPSG:" + crs, layer_name, "memory")
        line_provider = line_layer.dataProvider()

        fields = QgsFields()
        fields.append(QgsField("LINE_ID", QVariant.Int))

        for i in range(len(data_list[0]) - 3):
            fields.append(QgsField(f"data_{i}", QVariant.Double))

        line_provider.addAttributes(fields)
        line_layer.updateFields()

        line_features = {}

        for data in data_list:
            x, y, *values, line_id = data
            if line_id not in line_features:
                line_features[line_id] = []
            line_features[line_id].append((x, y, values))

        for line_id, points in line_features.items():
            coords = [QgsPointXY(x, y) for x, y, _ in points]
            geometry = QgsGeometry.fromPolylineXY(coords)

            feature = QgsFeature()
            feature.setGeometry(geometry)

            first_values = points[0][2]
            feature.setAttributes([line_id] + first_values)
            line_provider.addFeature(feature)

        QgsProject.instance().addMapLayer(line_layer)

        # Create the point layer
        point_layer = QgsVectorLayer(
            "Point?crs=EPSG:" + crs, f"{layer_name}_points", "memory"
        )
        point_provider = point_layer.dataProvider()

        fields = QgsFields()
        fields.append(QgsField("LINE_ID", QVariant.Int))

        for i in range(len(data_list[0]) - 3):
            fields.append(QgsField(f"data_{i}", QVariant.Double))

        point_provider.addAttributes(fields)
        point_layer.updateFields()

        for data in data_list:
            x, y, *values, line_id = data
            point = QgsGeometry.fromPointXY(QgsPointXY(x, y))

            feature = QgsFeature()
            feature.setGeometry(point)
            feature.setAttributes([line_id] + values)
            point_provider.addFeature(feature)

        if not point_layer.isValid():
            raise ValueError(f"Failed to load layer: {XYZ_file}")
        self.save_layer_as_shapefile(point_layer, output_path)

        layer2 = QgsVectorLayer(output_path, basename)
        # Add the layer to the current QGIS project
        QgsProject.instance().addMapLayer(layer2)

        layer_tree = QgsProject.instance().layerTreeRoot()
        layer_tree.findLayer(layer2.id()).setItemVisibilityChecked(False)

    def convert_RGB_to_grey(self, RGBGridPath, LUT):
        """
        Converts a 3-band RGB GeoTIFF to a grayscale GeoTIFF using a specified LUT (Look-Up Table).
        Args:
            RGBGridPath (str): The file path to the input 3-band RGB GeoTIFF.
            LUT (str): A comma-separated string of CSS color values representing the Look-Up Table.
        Returns:
            tuple:
                - result (bool): True if the conversion was successful, False otherwise.
                - RGBGridPath_gray (str or int): The file path to the generated grayscale GeoTIFF if successful,
                  or an error code (-3) if the operation failed.
        Notes:
            - The input RGB GeoTIFF must have at least 3 bands.
            - White ([255, 255, 255]) and black ([0, 0, 0]) pixels in the input are set to NaN in the output.
            - The LUT is parsed and used to map RGB values to scalar values via nearest neighbor lookup.
            - The output grayscale values are scaled based on user-defined minimum and maximum values.
            - The output GeoTIFF retains the geotransform and projection of the input dataset.
        Raises:
            Exception: If the output dataset cannot be created, possibly due to missing projection information.
        """

        result = False

        # Open the 3-band TIF using GDAL
        dataset = gdal.Open(RGBGridPath, gdal.GA_ReadOnly)
        if not dataset:
            self.iface.messageBar().pushMessage(
                "Unable to open the dataset.", level=Qgis.Warning, duration=15
            )
            return False, -3

        if dataset.RasterCount < 3:
            self.iface.messageBar().pushMessage(
                "Data file must have at least 3 layers", level=Qgis.Warning, duration=15
            )
            return False, -3

        red = dataset.GetRasterBand(1).ReadAsArray().astype(float)
        green = dataset.GetRasterBand(2).ReadAsArray().astype(float)
        blue = dataset.GetRasterBand(3).ReadAsArray().astype(float)

        transform = dataset.GetGeoTransform()
        projection = dataset.GetProjection()

        # Stack bands into an RGB array
        rgb_raster = np.dstack((red, green, blue))

        # Parse the LUT

        LUT = LUT.replace(" ", "")
        lut = self.parse_lut_string(LUT, num_entries=1024)

        """css_color_list = LUT.split(",")
        css_color_list.reverse()

        lut = self.generate_rgb_lut(css_color_list, num_entries=1024)"""
        if not lut:
            print("Couldn't generate LUT")
            return False, False

        scalar_values, lut_colors = zip(*lut)
        lut_colors = np.array(lut_colors) / 255.0  # Normalize LUT colors

        # Identify white and black pixels
        white_mask = (rgb_raster == [255, 255, 255]).all(axis=2)
        black_mask = (rgb_raster == [0, 0, 0]).all(axis=2)

        # Normalize raster RGB values to [0, 1]
        normalized_rgb = rgb_raster / 255.0

        # Flatten RGB raster for KDTree query
        reshaped_rgb = normalized_rgb.reshape(-1, 3)

        # Build a KDTree for nearest neighbor lookup
        lut_tree = cKDTree(lut_colors)
        distances, indices = lut_tree.query(reshaped_rgb)

        # Map nearest LUT color to scalar values
        scalar_grid = np.array(scalar_values)[indices].reshape(rgb_raster.shape[:2])

        # Set white and black areas to NaN
        scalar_grid[white_mask] = np.nan
        scalar_grid[black_mask] = np.nan

        # Scale data
        LUT_min = self.dlg.mQgsDoubleSpinBox_LUT_min.value()
        LUT_max = self.dlg.mQgsDoubleSpinBox_LUT_max.value()
        scalar_grid = (scalar_grid * (LUT_max - LUT_min)) + LUT_min

        # Prepare output file
        driver = gdal.GetDriverByName("GTiff")
        RGBGridPath_gray = self.insert_text_before_extension(RGBGridPath, "_gray")
        try:
            output_dataset = driver.Create(
                RGBGridPath_gray,
                dataset.RasterXSize,
                dataset.RasterYSize,
                1,
                gdal.GDT_Float32,
            )
        except:
            print(
                RGBGridPath_gray, dataset.RasterXSize, dataset.RasterYSize, projection
            )
            self.iface.messageBar().pushMessage(
                "Unable to create the output dataset, maybe check projection is set?",
                level=Qgis.Warning,
                duration=15,
            )
            return False, -3

        output_dataset.SetGeoTransform(transform)
        output_dataset.SetProjection(projection)

        # Write the scaled scalar grid to the output file
        output_band = output_dataset.GetRasterBand(1)
        output_band.WriteArray(scalar_grid)
        output_band.SetNoDataValue(np.nan)

        # Cleanup
        output_band = None
        output_dataset = None
        dataset = None

        result = True
        return result, RGBGridPath_gray

    def parse_lut_string(self, LUT, num_entries=1024):
        """
        Parse a LUT string and generate the appropriate RGB LUT.

        Parameters:
            LUT (str): String containing either CSS color names or RGB triplets
                    Examples:
                    - "red,green,blue"
                    - "(255,0,0),(0,255,0),(0,0,255)"
                    - "255,0,0,0,255,0,0,0,255"
                    - "(1.0,0,0),(0,1.0,0),(0,0,1.0)"
            num_entries (int): Number of entries in the LUT

        Returns:
            list: LUT list or False if parsing fails
        """
        import re

        # Remove spaces
        LUT = LUT.replace(" ", "")

        # Try to detect format and parse accordingly

        # Check for parentheses format: (R,G,B),(R,G,B)
        if "(" in LUT:
            rgb_pattern = r"\(([^)]+)\)"
            matches = re.findall(rgb_pattern, LUT)

            if matches:
                rgb_list = []
                for match in matches:
                    values = match.split(",")
                    try:
                        # Convert to floats (works for both int and float values)
                        if len(values) >= 3:
                            rgb = tuple(float(v) for v in values[:3])
                            rgb_list.append(rgb)
                        else:
                            raise ValueError("Not enough values for RGB")
                    except ValueError:
                        # Failed to parse as numbers, fall back to CSS colors
                        break
                else:
                    # Successfully parsed all as RGB
                    rgb_list.reverse()
                    return self.generate_rgb_lut_from_rgb(rgb_list, num_entries)

        # Try flat number format or CSS colors
        elements = LUT.split(",")

        # First, try to parse as numbers
        try:
            numbers = [float(elem) for elem in elements]

            # Check if we have RGB triplets (multiple of 3)
            if len(numbers) >= 3 and len(numbers) % 3 == 0:
                # Group into RGB triplets
                rgb_list = []
                for i in range(0, len(numbers), 3):
                    rgb = tuple(numbers[i : i + 3])
                    rgb_list.append(rgb)

                rgb_list.reverse()
                return self.generate_rgb_lut_from_rgb(rgb_list, num_entries)
        except ValueError:
            # Not all numbers, continue to CSS parsing
            pass

        # Default to CSS color names
        css_color_list = elements
        css_color_list.reverse()
        return self.generate_rgb_lut(css_color_list, num_entries)

    def generate_rgb_lut_from_rgb(self, rgb_list, num_entries=1024):
        """
        Generate an RGB LUT list from a list of RGB triplets.

        Parameters:
            rgb_list (list): List of RGB triplets. Each triplet can be:
                            - (R, G, B) where R, G, B are 0-255 integers
                            - (R, G, B) where R, G, B are 0-1 floats
                            - [R, G, B] lists work too
            num_entries (int): Total number of entries in the LUT.

        Returns:
            list: List of [decimal index, (R, G, B)] where R, G, B are 0-255 integers.
        """
        try:
            import matplotlib.colors as mcolors
        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for some functions, but not installed: matplotlib\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install matplotlib",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False

        # Normalize RGB values to 0-1 range if needed
        normalized_rgb_list = []
        for rgb in rgb_list:
            r, g, b = rgb[0], rgb[1], rgb[2]

            # Check if values are in 0-255 range (integers or floats > 1)
            if any(val > 1 for val in [r, g, b]):
                # Normalize from 0-255 to 0-1
                normalized_rgb = (r / 255.0, g / 255.0, b / 255.0)
            else:
                # Already in 0-1 range
                normalized_rgb = (r, g, b)

            normalized_rgb_list.append(normalized_rgb)

        # Normalize indices to decimal values between 0 and 1
        decimal_indices = np.linspace(0, 1, num_entries)

        # Create a continuous colormap using the normalized RGB list
        try:
            cmap = mcolors.LinearSegmentedColormap.from_list(
                "custom_cmap", normalized_rgb_list
            )
        except:
            return False

        # Generate RGB values for each index
        rgb_colors = [cmap(i)[:3] for i in decimal_indices]
        rgb_colors_255 = [
            (int(r * 255), int(g * 255), int(b * 255)) for r, g, b in rgb_colors
        ]

        # Combine decimal indices and RGB tuples
        lut = [
            [round(decimal_index, 6), rgb]
            for decimal_index, rgb in zip(decimal_indices, rgb_colors_255)
        ]

        return lut

    def generate_rgb_lut(self, css_color_list, num_entries=1024):
        """
        Generate an RGB LUT list from a list of CSS color names.

        Parameters:
            css_color_list (list): List of CSS color names recognized by Matplotlib.
            num_entries (int): Total number of entries in the LUT.

        Returns:
            list: List of [decimal index, (R, G, B)] where R, G, B are 0-255 integers.
        """
        try:
            import matplotlib.colors as mcolors

        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for some functions, but not installed: matplotlib\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install matplotlib",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False

        # Normalize indices to decimal values between 0 and 1
        decimal_indices = np.linspace(0, 1, num_entries)

        # Create a continuous colormap using the input CSS color list
        try:
            cmap = mcolors.LinearSegmentedColormap.from_list(
                "custom_cmap", css_color_list
            )
        except:

            return False

        # Generate RGB values for each index
        rgb_colors = [cmap(i)[:3] for i in decimal_indices]
        rgb_colors_255 = [
            (int(r * 255), int(g * 255), int(b * 255)) for r, g, b in rgb_colors
        ]

        # Combine decimal indices and RGB tuples
        lut = [
            [round(decimal_index, 6), rgb]
            for decimal_index, rgb in zip(decimal_indices, rgb_colors_255)
        ]

        return lut

    def vector_layer_to_dataframe(self, layer, attribute_name):
        """
        Converts a vector layer to a pandas DataFrame containing x, y coordinates and a chosen attribute.

        Parameters:
        - layer_name (str): The name of the vector layer in the QGIS project.
        - attribute_name (str): The attribute to include in the DataFrame.

        Returns:
        - pd.DataFrame: A DataFrame with columns 'x', 'y', and the chosen attribute.
        """
        # Get the layer by name

        # Check if the attribute exists in the layer
        if attribute_name not in [field.name() for field in layer.fields()]:
            raise ValueError(f"Attribute '{attribute_name}' not found in layer")

        # Extract features and build DataFrame
        data = self.extract_features_to_array(layer, attribute_name)

        crs = layer.crs()  # Get the CRS of the layer
        epsg_code = crs.postgisSrid()  # Retrieve the EPSG code

        # Create and return a DataFrame
        return data, epsg_code

    def extract_features_to_array(self, layer, attribute_name):
        """
        Extract x, y, value data from a QGIS layer into a 3×n NumPy array.

        :param layer: QGIS layer to extract features from.
        :param attribute_name: The attribute name to extract as the value.
        :return: 3×n NumPy array of x, y, and value data.
        """
        data = []

        for feature in layer.getFeatures():
            geom = feature.geometry()
            if geom and geom.isMultipart():
                # Handle multipart geometries
                for part in geom.asMultiPoint():
                    data.append([part.x(), part.y(), feature[attribute_name]])
            elif geom:
                # Handle single-part geometries
                point = geom.asPoint()
                data.append([point.x(), point.y(), feature[attribute_name]])

        # Convert to a NumPy array and return as a 3×n array
        return np.array(data)

    def create_temp_raster_mask_from_convex_hull(
        self, vector_layer, extent, cell_size=10
    ):
        """
        Creates a temporary raster mask based on the convex hull of a vector layer of points.

        :param vector_layer: The input vector layer (points).
        :param cell_size: Cell size for the output raster (default: 10).
        :return: Path to the temporary raster mask.
        """
        feedback = QgsProcessingFeedback()
        try:
            # Step 1: Generate convex hull
            convex_hull_output = processing.run(
                "qgis:minimumboundinggeometry",
                {"INPUT": vector_layer, "TYPE": 3, "OUTPUT": "TEMPORARY_OUTPUT"},
                feedback=feedback,
            )
            convex_hull_layer = convex_hull_output["OUTPUT"]

            # Step 2: Rasterize the convex hull into a temporary file
            temp_raster_path = tempfile.NamedTemporaryFile(
                suffix=".tif", delete=False
            ).name
            rasterize_output = processing.run(
                "gdal:rasterize",
                {
                    "INPUT": convex_hull_layer,
                    "FIELD": None,  # Use the entire polygon
                    "BURN": 1,  # Value to burn into the raster
                    "UNITS": 1,  # Cell size in map units
                    "WIDTH": cell_size,
                    "HEIGHT": cell_size,
                    "EXTENT": extent,
                    "NODATA": 0,  # Value for no data cells
                    "OUTPUT": temp_raster_path,
                },
                feedback=feedback,
            )

            print(f"Temporary raster mask created: {temp_raster_path}")
            return rasterize_output["OUTPUT"]

        except Exception as e:
            print(f"Error creating raster mask: {e}")
            return None
        # replace with specific processor calls so raster clipping can be done easily...

    def procWTMM(self):
        """
        Processes the WTMM (Wavelet Transform Modulus Maxima) for the selected layer in QGIS.
        This function performs the following steps:
        1. Checks for the required Python packages (`pywt` and `matplotlib`) and prompts the user to install them if missing.
        2. Retrieves the selected vector layer and validates its geometry type (line or point).
        3. Depending on the geometry type:
           - For line layers, retrieves data from the profile.
           - For point layers, regularizes selected points and calculates spacing.
        4. Computes the WTMM 1D analysis using the `WTMM` class.
        5. Visualizes the results, including the D(h) vs h spectrum and WTMM 1D visualization.
        Returns:
            list: The results of the WTMM 1D analysis if successful.
            None: If the selected layer is invalid or data retrieval fails.
            False: If required Python packages are missing.
        Raises:
            ImportError: If the required packages (`pywt` or `matplotlib`) are not installed.
        Notes:
            - This function is designed to work within the QGIS environment.
            - The user is expected to select appropriate layers and configure parameters via the plugin interface.
        """

        try:
            import pywt
        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for some functions, but not installed: PyWavelets\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install PyWavelets",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False

        try:
            import matplotlib.pyplot as plt

        except ImportError:
            QMessageBox.information(
                None,  # Parent widget
                "",
                "Missing Packages for SGTool: "  # Window title
                + f"The following Python packages are required for some functions, but not installed: matplotlib\n\n"
                "Please open the QGIS Python Console and run the following command:\n\n"
                f"!pip3 install matplotlib",  # Message text
                QMessageBox_Ok,  # Buttons parameter
            )
            return False

        line_layer_name = self.dlg.mMapLayerComboBox_selectVectors.currentText()
        data_name = self.dlg.mFieldComboBox_data.currentText()
        if line_layer_name.split(".")[-1].lower() == "xyz":
            data_name = "data_2"
        line_layer = QgsProject.instance().mapLayersByName(line_layer_name)[0]

        if line_layer is None:
            print("No layer selected")
            return None

        if (
            isinstance(line_layer, QgsVectorLayer)
            and line_layer.geometryType() == QgsWkbTypes.LineGeometry
            and self.dlg.mMapLayerComboBox_selectGrid_worms.currentText() != ""
        ):
            raster_layer_name = (
                self.dlg.mMapLayerComboBox_selectGrid_worms.currentText()
            )
            data = self.get_data_from_profile()
            if data is None:
                print("No data retrieved from profile")
                return None
            plot_layer_name = raster_layer_name
        elif (
            isinstance(line_layer, QgsVectorLayer)
            and line_layer.geometryType() == QgsWkbTypes.PointGeometry
        ):
            spacing = self.dlg.doubleSpinBox_wtmm_spacing.value()
            if self.dlg.doubleSpinBox_wtmm_spacing.value() == 0:
                spacing = "auto"

            new_coords, data, median_spacing = self.regularize_selected_points(
                line_layer.name(), data_name, spacing=spacing, num_points=None
            )
            plot_layer_name = line_layer_name

        else:
            print("Selected layer is not a valid line or point layer.")
            return None

        wtmm = WTMM()
        results = wtmm.wtmm_1d(
            data,
            wavelet="mexh",
            num_scales=15,
            threshold_rel=0.05,  # Lower threshold to detect more maxima
            min_distance=3,
        )
        if results:

            # Plot the D(h) vs h spectrum
            fig, ax = plt.subplots(figsize=(8, 6))

            wtmm.plot_Dh_vs_h(data, plot_layer_name, results, ax=ax)
            plt.tight_layout()
            plt.show()
            wtmm.visualize_wtmm_1d(
                data,
                plot_layer_name,
                results,
                line_number=int(self.dlg.mFieldComboBox_feature.currentText()),
                save_path=None,
            )

            plt.show()
            return results

    def get_data_from_profile(self):
        """
        Extracts data from a raster layer along a selected line feature in a QGIS project.
        This method performs the following steps:
        1. Retrieves the selected raster layer and line layer from the QGIS interface.
        2. Validates the input parameters, including line spacing and selected layers.
        3. Generates points along the selected line feature at specified intervals.
        4. Samples raster values at the generated points.
        5. Returns the sampled raster values as a NumPy array.
        Returns:
            numpy.ndarray: A NumPy array containing the sampled raster values along the line feature.
        Notes:
            - The line spacing must be non-zero and ideally greater than the grid cell size.
            - If no raster layer is explicitly selected, the method attempts to find one in the project.
            - The method uses QGIS's native "Points Along Lines" processing algorithm to generate points.
        Raises:
            None: The method handles errors internally and returns None if any issue occurs.
        """

        raster_layer_name = self.dlg.mMapLayerComboBox_selectGrid_worms.currentText()
        line_spacing = float(self.dlg.doubleSpinBox_wtmm_spacing.value())
        if line_spacing == 0:
            self.iface.messageBar().pushMessage(
                "Spacing must be non-zero for extraciton from grid, and ideally greater than grid cell size",
                level=Qgis.Warning,
                duration=15,
            )
            return
        line_layer_name = self.dlg.mMapLayerComboBox_selectVectors.currentText()
        line_layer = QgsProject.instance().mapLayersByName(line_layer_name)[0]
        line_layer.selectByIds([int(self.dlg.mFieldComboBox_feature.currentText())])

        if line_layer is None:
            print("No layer selected")
            return None

        if (
            isinstance(line_layer, QgsVectorLayer)
            and line_layer.geometryType() == QgsWkbTypes.LineGeometry
        ):

            # Find the specified raster layer
            raster_layer = None
            if raster_layer_name:
                layers = QgsProject.instance().mapLayersByName(raster_layer_name)
                if layers:
                    raster_layer = layers[0]
            else:
                # If no name provided, try to find a raster layer in the project
                for layer in QgsProject.instance().mapLayers().values():
                    if isinstance(layer, QgsRasterLayer):
                        raster_layer = layer
                        break

            if not raster_layer:
                print("No raster layer found")
                return None

            # Use the processing algorithm to create points along the geometry
            # We'll use a temporary layer for the points
            params = {
                "INPUT": QgsProcessingFeatureSourceDefinition(
                    line_layer.id(),
                    selectedFeaturesOnly=True,
                    featureLimit=-1,
                    geometryCheck=QgsFeatureRequest.GeometryAbortOnInvalid,
                ),
                "DISTANCE": line_spacing,
                "START_OFFSET": 0,
                "END_OFFSET": 0,
                "OUTPUT": "memory:",
            }

            # Run the algorithm
            result = processing.run("native:pointsalonglines", params)
            points_layer = result["OUTPUT"]

            # Get the total number of points
            num_points = points_layer.featureCount()

            # Create arrays to store the points and raster values
            points_array = np.zeros((num_points, 2))
            raster_values = np.zeros(num_points)
            distances = np.zeros(num_points)

            # Get the raster data provider
            provider = raster_layer.dataProvider()

            # Loop through all points and sample the raster
            for i, point_feature in enumerate(points_layer.getFeatures()):
                # Get the point geometry
                point_geom = point_feature.geometry()
                point = point_geom.asPoint()

                # Get the distance attribute (added by the processing algorithm)
                distance = point_feature["distance"]

                # Add to our arrays
                points_array[i] = [point.x(), point.y()]

                # Sample the raster at this point (band 1)
                value, valid = provider.sample(QgsPointXY(point.x(), point.y()), 1)

                # Store the raster value and distance
                raster_values[i] = value if valid else np.nan
                distances[i] = distance

            # Clean up temporary layer
            QgsProject.instance().removeMapLayer(points_layer.id())

            # Create a complete result with distances, coordinates, and raster values
            result = {
                "distances": distances,
                "points": points_array,
                "values": raster_values,
            }
            data = np.array(result["values"])

            return data

    def regularize_selected_points(
        self, layer_name, field_name, spacing="auto", num_points=None
    ):
        """
        Convert selected points from a QGIS point layer to a NumPy array
        and regularize their spacing along a line

        Parameters:
        layer_name (str): Name of the point layer in the QGIS project
        field_name (str): Name of the field to extract
        spacing (float or "auto", optional): Desired spacing between points or "auto" to use median spacing
        num_points (int, optional): Desired number of points. If None and spacing is None,
                                original number of points is used

        Returns:
        tuple: (numpy.ndarray of coordinates, numpy.ndarray of field values, float: spacing used)
        """
        # Get the layer by name
        layer = QgsProject.instance().mapLayersByName(layer_name)[0]
        selection_value = self.dlg.mFieldComboBox_feature.currentText()
        if "LINE_ID" in [field.name() for field in layer.fields()]:
            layer.selectByExpression("LINE_ID = " + selection_value)
        else:
            layer.selectByExpression("LINE = " + selection_value)

        if not layer:
            print(f"Layer '{layer_name}' not found")
            return None, None, None

        # Check if the field exists in the layer
        fields = layer.fields()
        field_idx = fields.indexFromName(field_name)

        if field_idx == -1:
            print(f"Field '{field_name}' not found in layer '{layer_name}'")
            return None, None, None

        # Get the selected features
        selected_features = layer.selectedFeatures()

        if not selected_features:
            print(f"No features selected in layer '{layer_name}'")
            return None, None, None

        # Get the coordinates and field values from selected features
        coords = []
        field_values = []

        for feature in selected_features:
            geom = feature.geometry()
            point = geom.asPoint()
            coords.append((point.x(), point.y()))
            field_values.append(feature[field_name])

        # Convert to numpy arrays
        coords = np.array(coords)
        field_values = np.array(field_values)

        # Determine the dominant axis (x or y) by checking the range
        x_range = np.max(coords[:, 0]) - np.min(coords[:, 0])
        y_range = np.max(coords[:, 1]) - np.min(coords[:, 1])

        # Sort the points based on the dominant axis
        if x_range >= y_range:
            # Sort by x-coordinate
            print("Sorting points along X axis (dominant)")
            sort_idx = np.argsort(coords[:, 0])
        else:
            # Sort by y-coordinate
            print("Sorting points along Y axis (dominant)")
            sort_idx = np.argsort(coords[:, 1])

        coords_sorted = coords[sort_idx]
        field_values_sorted = field_values[sort_idx]

        # Calculate distances between consecutive points
        distances = [0]
        point_distances = []

        for i in range(1, len(coords_sorted)):
            prev = coords_sorted[i - 1]
            curr = coords_sorted[i]
            d = np.sqrt((curr[0] - prev[0]) ** 2 + (curr[1] - prev[1]) ** 2)
            point_distances.append(d)
            distances.append(distances[-1] + d)

        total_length = distances[-1]

        # Determine regularization parameters
        if spacing == "auto":
            # Calculate median point spacing
            median_spacing = np.median(point_distances)
            print(f"Using auto-calculated median spacing: {median_spacing}")
            num_new_points = int(total_length / median_spacing) + 1
            used_spacing = median_spacing
        elif isinstance(spacing, (int, float)):
            num_new_points = int(total_length / spacing) + 1
            used_spacing = spacing
        elif num_points is not None:
            num_new_points = num_points
            used_spacing = (
                total_length / (num_new_points - 1) if num_new_points > 1 else 0
            )
        else:
            num_new_points = len(coords_sorted)
            used_spacing = (
                total_length / (num_new_points - 1) if num_new_points > 1 else 0
            )

        # Create evenly spaced points along the path
        new_distances = np.linspace(0, total_length, num_new_points)

        # Interpolate coordinates along the path
        x_interp = interp1d(distances, coords_sorted[:, 0], kind="linear")
        y_interp = interp1d(distances, coords_sorted[:, 1], kind="linear")

        # Interpolate field values along the path
        field_interp = interp1d(distances, field_values_sorted, kind="linear")

        # Calculate new coordinates and field values
        new_coords = np.column_stack((x_interp(new_distances), y_interp(new_distances)))
        new_field_values = field_interp(new_distances)

        return new_coords, new_field_values, used_spacing

    def create_regularized_layer(
        self,
        new_coords,
        new_field_values,
        field_name,
        output_layer_name="RegularizedPoints",
    ):
        """
        Create a new point layer with regularized points

        Parameters:
        new_coords (numpy.ndarray): Array of coordinates (x, y)
        new_field_values (numpy.ndarray): Array of field values
        field_name (str): Name of the field
        output_layer_name (str): Name for the output layer

        Returns:
        QgsVectorLayer: The created layer
        """
        # Create a new memory layer
        layer = QgsVectorLayer(
            f"Point?crs=EPSG:4326&field={field_name}:double",
            output_layer_name,
            "memory",
        )

        # Get ready to add features
        provider = layer.dataProvider()

        # Add features
        features = []
        for i in range(len(new_coords)):
            feat = QgsFeature()
            point = QgsPointXY(new_coords[i][0], new_coords[i][1])
            feat.setGeometry(QgsGeometry.fromPointXY(point))
            feat.setAttributes([float(new_field_values[i])])
            features.append(feat)

        provider.addFeatures(features)

        # Add layer to the project
        QgsProject.instance().addMapLayer(layer)

        return layer

    def selected_points_to_numpy(self, layer_name, field_name):
        """
        Convert selected points from a QGIS point layer to a NumPy array for a specific field

        Parameters:
        layer_name (str): Name of the point layer in the QGIS project
        field_name (str): Name of the field to extract

        Returns:
        tuple: (numpy.ndarray of coordinates, numpy.ndarray of field values)
        """
        # Get the layer by name
        layer = QgsProject.instance().mapLayersByName(layer_name)[0]

        if not layer:
            print(f"Layer '{layer_name}' not found")
            return None, None

        # Check if the field exists in the layer
        fields = layer.fields()
        field_idx = fields.indexFromName(field_name)

        if field_idx == -1:
            print(f"Field '{field_name}' not found in layer '{layer_name}'")
            return None, None

        # Get the selected features
        selected_features = layer.selectedFeatures()

        if not selected_features:
            print(f"No features selected in layer '{layer_name}'")
            return None, None

        # Get the coordinates and field values from selected features
        coords = []
        field_values = []

        for feature in selected_features:
            geom = feature.geometry()
            point = geom.asPoint()
            coords.append((point.x(), point.y()))
            field_values.append(feature[field_name])

        # Convert to numpy arrays
        coords = np.array(coords)
        field_values = np.array(field_values)

        return coords, field_values

    def safe_numeric_sort_key(self, value):
        """
        A safe sorting key function that can handle both integers and floats.
        This also handles cases where the value might be a string representation of a number.

        Args:
            value: The value to convert to a numeric type for sorting

        Returns:
            float or original value: The numeric value for sorting or the original value if conversion fails
        """
        if value is None:
            return float("-inf")  # None values come first in sorting

        # Try to convert to float first
        try:
            return float(value)
        except (ValueError, TypeError):
            # If it can't be converted to a number, return the original value
            return value

    def update_wavelet_choices(self):
        """
        Updates the wavelet choices in the user interface based on the selected vector layer.
        This method performs the following steps:
        1. Checks if a vector layer is selected in the `mMapLayerComboBox_selectVectors` combo box.
        2. Retrieves the selected vector layer and enables the `mFieldComboBox_data` combo box.
        3. Verifies if the selected layer contains a "LINE_ID" field. If not, the function exits.
        4. Depending on the geometry type of the layer:
            - For PointGeometry:
                - Clears the `mFieldComboBox_feature` combo box.
                - Collects unique "LINE_ID" values from the features in the layer.
                - Sorts the unique values and populates the `mFieldComboBox_feature` combo box.
                - Sets the layer for the `mFieldComboBox_data` combo box.
            - For LineGeometry:
                - Clears the `mFieldComboBox_feature` combo box.
                - Disables the `mFieldComboBox_data` combo box.
                - Collects unique feature IDs from the layer.
                - Sorts the unique IDs and populates the `mFieldComboBox_feature` combo box.
        Returns:
            None
        """

        if self.dlg.mMapLayerComboBox_selectVectors.currentText() != "":
            line_layer_name = self.dlg.mMapLayerComboBox_selectVectors.currentText()
            line_layer = QgsProject.instance().mapLayersByName(line_layer_name)
            if not line_layer:
                return

            line_layer = line_layer[0]

            self.dlg.mFieldComboBox_data.setEnabled(True)

            if line_layer.geometryType() == QgsWkbTypes.PointGeometry:
                # Check if the layer has the LINE_ID field
                if "LINE_ID" not in [
                    field.name() for field in line_layer.fields()
                ] and "LINE" not in [field.name() for field in line_layer.fields()]:
                    # Field doesn't exist, break out or handle the error
                    # QMessageBox.warning(None, "Field Missing", "The layer does not contain a LINE_ID field.")
                    return  # This will exit the current function

                self.dlg.mFieldComboBox_feature.clear()
                unique_values = []

                for feature in line_layer.getFeatures():
                    if "LINE_ID" in [field.name() for field in line_layer.fields()]:
                        value = str(feature["LINE_ID"])
                    else:
                        value = str(feature["LINE"])

                    if value not in unique_values:
                        unique_values.append(value)

                # Sort the values (optional)
                unique_values = sorted(unique_values, key=self.safe_numeric_sort_key)

                # Add to combo box
                self.dlg.mFieldComboBox_feature.addItems(unique_values)
                self.dlg.mFieldComboBox_data.setLayer(line_layer)
            elif line_layer.geometryType() == QgsWkbTypes.LineGeometry:
                self.dlg.mFieldComboBox_feature.clear()
                self.dlg.mFieldComboBox_data.setEnabled(False)
                unique_values = []

                for feature in line_layer.getFeatures():
                    value = feature.id()
                    if value not in unique_values:
                        unique_values.append(str(value))

                # Sort the values (optional)
                unique_values = sorted(unique_values, key=int)

                # Add to combo box
                self.dlg.mFieldComboBox_feature.addItems(unique_values)


class OddPositiveIntegerValidator(QValidator):
    def validate(self, input_text, pos):
        """
        Validates the text input to allow only odd positive integers.
        """
        if not input_text:  # Allow empty input (to clear the field)
            return QValidator.Intermediate, input_text, pos

        if not input_text.isdigit():  # Only digits are allowed
            return QValidator.Invalid, input_text, pos

        value = int(input_text)
        if value > 0 and value % 2 == 1:  # Check for positive odd numbers
            return QValidator.State.Acceptable, input_text, pos
        else:
            return QValidator.State.Intermediate, input_text, pos

    def fixup(self, input_text):
        """
        Corrects invalid input to the nearest odd positive integer.
        """
        try:
            value = int(input_text)
            if value <= 0:  # Make it a positive number
                return "1"
            elif value % 2 == 0:  # Make it odd
                return str(value + 1)
            else:
                return input_text
        except ValueError:
            return "1"  # Default to 1 if the input cannot be converted
