"""
/***************************************************************************
 qgSurf - plugin for Quantum GIS

 Processing of geological planes and surfaces

                              -------------------
        begin                : 2011-12-21
        copyright            : (C) 2011-2026 by Mauro Alberti
        email                : alberti.m65@gmail.com

 ***************************************************************************/

# licensed under the terms of GNU GPL 3

/***************************************************************************
 *                                                                         *
 *   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 3 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""

from math import atan, isfinite
import math

from qgis.PyQt.QtWidgets import *
from qgis.PyQt.QtCore import QSignalBlocker
from qgis.core import QgsMessageLog
from qgis.core import Qgis

try:
    from osgeo import ogr
except ImportError:
    import ogr

from geogst.core.geometries.grids.rasters import *
from geogst.io.rasters.gdal import *

from ..apps.qgis.canvas import *
from ..apps.qgis.orientations import *
from ..apps.qgis.points import *
from ..apps.qgis.project import *
from ..apps.qgis.projections import *
from ..apps.qt.messages import *

from ..configurations.colors import *
from ..configurations.fields import *
from ..configurations.general_params import *
from ..tools.config.constants import (
    BufferConstants,
    UIConstants,
    GeometryConstants,
    EPSGCodes,
)


class DemPlaneIntersectionWidget(QWidget):
    """
    Constructor
    
    """

    line_colors = [
        "white",
        "red",
        "blue",
        "yellow",
        "orange",
        "brown"
    ]

    required_text = '--  required  --'

    def __init__(self,
                 plugin_folder: str,
                 module_name: str,
                 canvas,
                 plugin_qaction):

        super(DemPlaneIntersectionWidget, self).__init__()

        self.plugin_folder = plugin_folder
        self.module_name = module_name

        self.canvas = canvas
        self.plugin = plugin_qaction

        self.setup_gui()
        self.initialize_parameters()
        self.update_crs_settings()

    def initialize_parameters(self):

        self.previous_tool = None
        self.use_elevation_from_dem = False
        self.map_tool = None

        self.reset_everything()

    def reset_post_dem_change(self):

        self.disable_canvas_tools()
        self.reset_input_dem_parameters()
        self.reset_intersections()

    def reset_everything(self):

        self.disable_canvas_tools()
        self.reset_all_io_parameters()

    def reset_source_point_parameters(self):

        self.reset_source_point_in_map_QPushButton.setEnabled(False)
        self.reset_source_point_and_intersections()

    def reset_all_io_parameters(self):

        self.reset_all_input_parameters()
        self.reset_intersections()

    def reset_all_input_parameters(self):

        self.reset_input_dem_parameters()
        self.reset_source_point()
        self.reset_input_geoplane_parameters()
        self.reset_buffer_parameters()

    def reset_buffer_parameters(self):
        """Reset buffer zone parameters to default values."""
        self.use_buffer_zone_QCheckBox.setChecked(False)
        self.buffer_radius_QDoubleSpinBox.setValue(BufferConstants.DEFAULT_BUFFER_RADIUS_METERS)
        self.buffer_radius_QDoubleSpinBox.setEnabled(False)

    def reset_input_dem_parameters(self):

        self.dem = None
        self.grid = None
        self.current_source_point_elevation_value = None

    def disable_canvas_tools(self):

        self.set_source_point_in_map_QPushButton.setEnabled(False)
        self.reset_source_point_in_map_QPushButton.setEnabled(False)

        try:
            self.map_tool.canvasClicked.disconnect(self.update_source_point_location)
        except:
            pass

        try:
            self.restore_previous_map_tool(self.map_tool)
        except:
            pass

    def reset_source_point_and_intersections(self):

        self.reset_intersections()
        self.reset_source_point()

    def reset_intersections(self):

        if hasattr(self, 'intersections_markers_list'):
            self.remove_intersection_markers_from_canvas()
        self.reset_intersections_markers_and_values()

    def reset_intersections_markers_list(self):

        self.intersections_markers_list = []

    def reset_intersections_values(self):

        self.intersections_x = []
        self.intersections_y = []
        self.intersections_xprt = {}
        self.inters = None

    def reset_intersections_markers_and_values(self):

        self.reset_intersections_markers_list()
        self.reset_intersections_values()

    def reset_source_point(self):

        self.clear_source_point_qlineedit_forms()
        if hasattr(self, 'source_point_marker_list'):
            self.remove_source_point_marker_from_canvas()
        self.reset_source_point_coordinates()
        self.reset_source_point_marker_list()

    def reset_source_point_coordinates(self):

        self.source_point_longitude = None
        self.source_point_latitude = None
        self.source_point_elevation = None

    def reset_source_point_marker_list(self):

        self.source_point_marker_list = []

    '''
    def remove_intersection_markers_from_canvas(self):

        for mrk in self.intersections_markers_list:
            self.canvas.scene().removeItem(mrk)

    def remove_source_point_marker_from_canvas(self):

        for mrk in self.source_point_marker_list:
            self.canvas.scene().removeItem(mrk)
    '''

    def _safe_remove_marker(self, mrk):
        """
        Rimuove il marker dalla sua scena solo se esiste ancora e appartiene a una scena.
        Evita i warning QGraphicsScene::removeItem.
        """
        try:
            if mrk is None:
                return

            scene = mrk.scene()
            if scene is not None:
                scene.removeItem(mrk)
            # opzionale: nascondi comunque il marker
            mrk.setVisible(False)

        except Exception:
            # se il wrapped C++ object è già distrutto o altro, ignora tranquillamente
            pass

    def remove_intersection_markers_from_canvas(self):

        if not hasattr(self, 'intersections_markers_list'):
            return

        for mrk in self.intersections_markers_list:
            self._safe_remove_marker(mrk)

        # importantissimo: svuotare la lista
        self.intersections_markers_list = []

    def remove_source_point_marker_from_canvas(self):

        if not hasattr(self, 'source_point_marker_list'):
            return

        for mrk in self.source_point_marker_list:
            self._safe_remove_marker(mrk)

        # idem
        self.source_point_marker_list = []

    def reset_input_geoplane_parameters(self):

        self.source_dip_direction = None
        self.source_dip_angle = None

    def refresh_raster_layer_list(self):

        # Use QSignalBlocker to prevent signal emissions during combobox update
        with QSignalBlocker(self.choose_source_dem_QComboBox):
            try:
                self.choose_source_dem_QComboBox.currentIndexChanged[int].disconnect(self.get_chosen_dem)
            except:
                pass

            try:
                self.reset_everything()
            except:
                pass

            try:
                self.choose_source_dem_QComboBox.clear()
            except:
                return

            self.choose_source_dem_QComboBox.addItem(self.required_text)

            self.raster_layers = loaded_raster_layers()
            if self.raster_layers is None or len(self.raster_layers) == 0:
                return
            for layer in self.raster_layers:
                self.choose_source_dem_QComboBox.addItem(layer.name())

            self.choose_source_dem_QComboBox.currentIndexChanged[int].connect(self.get_chosen_dem)

    def setup_gui(self):

        layout = QVBoxLayout()

        widget = QTabWidget()

        widget.addTab(self.setup_point_intersections_tab(), "Processing")
        widget.addTab(self.setup_help_tab(), "Help")
        
        layout.addWidget(widget)

        self.setLayout(layout)
        self.adjustSize()

        self.setWindowTitle(f'{plugin_name} - {self.module_name}')

    def setup_point_intersections_tab(self):
        
        widget = QWidget()

        layout = QVBoxLayout()
        layout.addWidget(self.setup_source_dem())
        layout.addWidget(self.setup_tabs())

        widget.setLayout(layout)
        return widget

    def setup_source_dem(self):

        widget = QWidget()

        layout = QGridLayout()

        self.refresh_raster_layer_list_QPushButton = QPushButton("Update raster layer list")
        self.refresh_raster_layer_list_QPushButton.clicked[bool].connect(self.refresh_raster_layer_list)
        self.refresh_raster_layer_list_QPushButton.setEnabled(True)
        layout.addWidget(self.refresh_raster_layer_list_QPushButton, 0, 0, 1, 3)

        layout.addWidget(QLabel("Choose current DEM layer"), 1, 0, 1, 2)
        self.choose_source_dem_QComboBox = QComboBox()
        self.choose_source_dem_QComboBox.addItem(self.required_text)
        layout.addWidget(self.choose_source_dem_QComboBox, 1, 2, 1, 1)

        self.refresh_raster_layer_list()

        # load geographical and geological parameters

        layout.addWidget(QLabel(self.tr("Load geographical and geological parameters from file:")), 2, 0, 1, 2)

        input_file_browser = QPushButton(".....")
        input_file_browser.clicked.connect(self.load_input_processing_parameters)
        layout.addWidget(input_file_browser, 2, 2, 1, 1)

        widget.setLayout(layout)
        return widget

    def setup_tabs(self):  

        widget = QWidget()

        layout = QVBoxLayout()

        toolbox = QToolBox()
        toolbox.addItem(self.setup_geographic_parameters_tab(), 'Geographic parameters')
        toolbox.addItem(self.setup_geologic_parameters_tab(), 'Geological parameters')
        toolbox.addItem(self.setup_calculations_tab(), 'Calculate')
        toolbox.addItem(self.setup_output_tab(), 'Output')

        layout.addWidget(toolbox)
        widget.setLayout(layout)

        return widget

    def setup_help_tab(self):
        
        widget = QWidget()
        layout = QVBoxLayout()

        # About section

        text_browser = QTextBrowser(widget)

        url_path = f"file:///{self.plugin_folder}/help/help_di.html"
        text_browser.setSource(QUrl(url_path))
        text_browser.setSearchPaths([f'{self.plugin_folder}/help/images/di'])

        layout.addWidget(text_browser)

        widget.setLayout(layout)

        return widget

    def setup_geographic_parameters_tab(self):
        
        widget = QWidget()

        layout = QGridLayout()

        self.set_source_point_in_map_QPushButton = QPushButton("Set source point in map")
        self.set_source_point_in_map_QPushButton.clicked[bool].connect(self.set_source_point_in_map)
        layout.addWidget(self.set_source_point_in_map_QPushButton, 0, 0, 1, 3)

        self.lock_z_value_to_dem_surface_QCheckBox = QCheckBox("Get elevation from DEM surface")
        self.lock_z_value_to_dem_surface_QCheckBox.setChecked(True)
        self.lock_z_value_to_dem_surface_QCheckBox.stateChanged[int].connect(self.update_source_point_elevation_value)
        layout.addWidget(self.lock_z_value_to_dem_surface_QCheckBox, 1, 0, 1, 3)

        layout.addWidget(QLabel("Longitude (EPSG: 4326)"), 2, 0, 1, 1)
        self.source_point_longitude_QLineEdit = QLineEdit()
        self.source_point_longitude_QLineEdit.textEdited.connect(self.update_source_point_location)
        layout.addWidget(self.source_point_longitude_QLineEdit, 2, 1, 1, 2)
              
        layout.addWidget(QLabel("Latitude (EPSG: 4326)"), 3, 0, 1, 1)
        self.source_point_latitude_QLineEdit = QLineEdit()
        self.source_point_latitude_QLineEdit.textEdited.connect(self.update_source_point_location)
        layout.addWidget(self.source_point_latitude_QLineEdit, 3, 1, 1, 2)
                      
        layout.addWidget(QLabel("Elevation"), 4, 0, 1, 1)
        self.source_point_elevation_QLineEdit = QLineEdit()
        self.source_point_elevation_QLineEdit.textEdited.connect(self.check_elevation_congruence_with_dem)
        layout.addWidget(self.source_point_elevation_QLineEdit, 4, 1, 1, 2)

        self.reset_source_point_in_map_QPushButton = QPushButton("Reset source point")
        self.reset_source_point_in_map_QPushButton.clicked[bool].connect(self.reset_source_point_parameters)
        layout.addWidget(self.reset_source_point_in_map_QPushButton, 5, 0, 1, 3)

        # Add buffer zone controls
        self.use_buffer_zone_QCheckBox = QCheckBox("Limit calculation to buffer zone around source point")
        self.use_buffer_zone_QCheckBox.setChecked(False)
        layout.addWidget(self.use_buffer_zone_QCheckBox, 6, 0, 1, 3)

        layout.addWidget(QLabel("Buffer radius (map units)"), 7, 0, 1, 1)
        self.buffer_radius_QDoubleSpinBox = QDoubleSpinBox()
        self.buffer_radius_QDoubleSpinBox.setRange(
            BufferConstants.BUFFER_RADIUS_MIN_METERS,
            BufferConstants.BUFFER_RADIUS_MAX_METERS
        )
        self.buffer_radius_QDoubleSpinBox.setDecimals(1)
        self.buffer_radius_QDoubleSpinBox.setSingleStep(BufferConstants.BUFFER_RADIUS_SINGLE_STEP_METERS)
        self.buffer_radius_QDoubleSpinBox.setValue(BufferConstants.DEFAULT_BUFFER_RADIUS_METERS)
        self.buffer_radius_QDoubleSpinBox.setEnabled(False)
        layout.addWidget(self.buffer_radius_QDoubleSpinBox, 7, 1, 1, 2)

        # Connect checkbox to enable/disable spinbox
        self.use_buffer_zone_QCheckBox.stateChanged[int].connect(self.toggle_buffer_controls)

        widget.setLayout(layout)
        
        return widget

    def update_crs_settings(self):

        self.map_crs = self.canvas.mapSettings().destinationCrs()

    def check_elevation_congruence_with_dem(self):
        
        if self.use_elevation_from_dem and float(self.source_point_elevation_QLineEdit.text()) != self.current_source_point_elevation_value:
            self.use_elevation_from_dem = False
            self.lock_z_value_to_dem_surface_QCheckBox.setChecked(False)
            
        self.current_source_point_elevation_value = float(self.source_point_elevation_QLineEdit.text())

    def toggle_buffer_controls(self):
        """Enable/disable buffer radius control based on checkbox state."""
        is_enabled = self.use_buffer_zone_QCheckBox.isChecked()
        self.buffer_radius_QDoubleSpinBox.setEnabled(is_enabled)

    def setup_geologic_parameters_tab(self):

        widget = QWidget()
        layout = QGridLayout()

        layout.addWidget(QLabel("Input dip azimuth is relative to:"), 0, 0, 1, 3)

        self.direction_azimuth_is_relative_to_geographic_north_QRadioButton = QRadioButton("Geographic North")
        self.direction_azimuth_is_relative_to_geographic_north_QRadioButton.setChecked(True)
        layout.addWidget(self.direction_azimuth_is_relative_to_geographic_north_QRadioButton, 1, 1, 1, 1)

        self.direction_azimuth_is_relative_to_map_top_QRadioButton = QRadioButton("Current map top (Y axis)")
        layout.addWidget(self.direction_azimuth_is_relative_to_map_top_QRadioButton, 1, 2, 1, 1)

        dip_dir_label = QLabel("Dip direction")
        dip_dir_label.setAlignment (Qt.AlignmentFlag.AlignCenter)       
        layout.addWidget(dip_dir_label, 2, 0, 1, 2)

        dip_ang_label = QLabel("Dip angle")
        dip_ang_label.setAlignment(Qt.AlignmentFlag.AlignCenter)       
        layout.addWidget(dip_ang_label, 2, 2, 1, 1)
        
        self.dip_direction_QDial = QDial()
        self.dip_direction_QDial.setRange(UIConstants.DIP_DIRECTION_MIN, UIConstants.DIP_DIRECTION_MAX)
        self.dip_direction_QDial.setPageStep(1)
        self.dip_direction_QDial.setProperty("value", UIConstants.DIP_DIRECTION_DEFAULT)
        self.dip_direction_QDial.setSliderPosition(UIConstants.DIP_DIRECTION_DEFAULT)
        self.dip_direction_QDial.setTracking(True)
        self.dip_direction_QDial.setOrientation(Qt.Orientation.Vertical)
        self.dip_direction_QDial.setWrapping(True)
        self.dip_direction_QDial.setNotchTarget(UIConstants.DIP_DIRECTION_NOTCH_TARGET)
        self.dip_direction_QDial.setNotchesVisible(True)
        self.dip_direction_QDial.valueChanged[int].connect(self.update_dip_direction_spinbox)
        layout.addWidget(self.dip_direction_QDial, 3, 0, 1, 2)
                
        self.DAngle_verticalSlider = QSlider()
        self.DAngle_verticalSlider.setRange(UIConstants.DIP_ANGLE_MIN, UIConstants.DIP_ANGLE_MAX)
        self.DAngle_verticalSlider.setProperty("value", UIConstants.DIP_ANGLE_DEFAULT)
        self.DAngle_verticalSlider.setOrientation(Qt.Orientation.Vertical)
        self.DAngle_verticalSlider.setInvertedAppearance(True)
        self.DAngle_verticalSlider.setTickPosition(QSlider.TickPosition.TicksBelow)
        self.DAngle_verticalSlider.setTickInterval(UIConstants.DIP_ANGLE_TICK_INTERVAL)
        self.DAngle_verticalSlider.valueChanged[int].connect(self.update_dip_angle_spinbox)
        layout.addWidget(self.DAngle_verticalSlider, 3, 2, 1, 1)

        self.dip_direction_QDoubleSpinBox = QDoubleSpinBox()
        self.dip_direction_QDoubleSpinBox.setRange(
            UIConstants.DIP_DIRECTION_SPINBOX_MIN,
            UIConstants.DIP_DIRECTION_SPINBOX_MAX
        )
        self.dip_direction_QDoubleSpinBox.setDecimals(3)
        self.dip_direction_QDoubleSpinBox.setSingleStep(UIConstants.DIP_DIRECTION_SPINBOX_SINGLE_STEP)
        self.dip_direction_QDoubleSpinBox.valueChanged[float].connect(self.update_dip_direction_slider)
        layout.addWidget(self.dip_direction_QDoubleSpinBox, 4, 0, 1, 2)
         
        self.dip_angle_QDoubleSpinBox = QDoubleSpinBox()
        self.dip_angle_QDoubleSpinBox.setRange(
            UIConstants.DIP_ANGLE_SPINBOX_MIN,
            UIConstants.DIP_ANGLE_SPINBOX_MAX
        )
        self.dip_angle_QDoubleSpinBox.setDecimals(3)
        self.dip_angle_QDoubleSpinBox.setSingleStep(UIConstants.DIP_ANGLE_SPINBOX_SINGLE_STEP)
        self.dip_angle_QDoubleSpinBox.setProperty("value", UIConstants.DIP_ANGLE_SPINBOX_DEFAULT)
        self.dip_angle_QDoubleSpinBox.valueChanged[float].connect(self.update_dip_angle_slider)
        layout.addWidget(self.dip_angle_QDoubleSpinBox, 4, 2, 1, 1)

        widget.setLayout(layout)

        return widget
        
    def setup_calculations_tab(self):

        widget = QWidget()
        layout = QVBoxLayout()

        dem_groupbox = QGroupBox(self.tr("Intersections of DEM with geological plane -> intersection points"))

        dem_layout = QGridLayout()

        dem_layout.addWidget(QLabel("Marker color:"), 1, 0, 1, 1)

        self.intersections_color_QComboBox = QComboBox()
        self.intersections_color_QComboBox.insertItems(
            0,
            INTERSECTIONS_COLOR
        )

        self.intersections_color_QComboBox.currentIndexChanged[int].connect(self.plot_intersections)

        dem_layout.addWidget(self.intersections_color_QComboBox, 1, 1, 1, 1)

        self.calculate_intersections_QPushButton = QPushButton("Calculate")
        self.calculate_intersections_QPushButton.clicked[bool].connect(self.calculate_intersections)
        dem_layout.addWidget(self.calculate_intersections_QPushButton, 2, 0, 1, 4)

        self.reset_intersections_QPushButton = QPushButton("Reset")
        self.reset_intersections_QPushButton.clicked.connect(self.reset_intersections)
        dem_layout.addWidget(self.reset_intersections_QPushButton, 3, 0, 1, 4)

        dem_groupbox.setLayout(dem_layout)

        layout.addWidget(dem_groupbox)

        # final setup

        widget.setLayout(layout)

        return widget

    def setup_output_tab(self):

        widget = QWidget()

        layout = QGridLayout()

        # intersections

        self.export_last_intersection_result_QPushButton = QPushButton("Save intersections as shapefile:")
        self.export_last_intersection_result_QPushButton.clicked.connect(self.export_last_intersection_result)
        layout.addWidget(self.export_last_intersection_result_QPushButton, 0, 0, 1, 2)

        self.output_shapefile_path_QLineEdit = QLineEdit()
        layout.addWidget(self.output_shapefile_path_QLineEdit, 0, 2, 1, 1)

        self.define_output_shapefile_path_QPushButton = QPushButton(".....")
        self.define_output_shapefile_path_QPushButton.clicked.connect(self.define_output_vector_file_path)
        layout.addWidget(self.define_output_shapefile_path_QPushButton, 0, 3, 1, 1)

        layout.addWidget(QLabel(self.tr("with geometry:")), 1, 0, 1, 1)

        save_qbuttongroup = QButtonGroup(widget)

        self.save_intersections_as_points_QRadioButton = QRadioButton("points")
        self.save_intersections_as_points_QRadioButton.setChecked(True)
        save_qbuttongroup.addButton(self.save_intersections_as_points_QRadioButton, 0)
        layout.addWidget(self.save_intersections_as_points_QRadioButton, 1, 1, 1, 1)

        self.save_intersections_as_lines_QRadioButton = QRadioButton("lines")
        self.save_intersections_as_lines_QRadioButton.setChecked(True)
        save_qbuttongroup.addButton(self.save_intersections_as_lines_QRadioButton, 0)
        layout.addWidget(self.save_intersections_as_lines_QRadioButton, 1, 2, 1, 1)

        self.load_output_in_map_QCheckBox = QCheckBox("load intersection shapefile in project")
        self.load_output_in_map_QCheckBox.setChecked(True)
        layout.addWidget(self.load_output_in_map_QCheckBox, 2, 0, 1, 2)

        # separation line

        line = QFrame()
        line.setFrameShape(QFrame.Shape.HLine)
        line.setFixedHeight(3)
        line.setFrameShadow(QFrame.Shadow.Raised)
        line.setSizePolicy(QSizePolicy.Policy.Expanding, QSizePolicy.Policy.Fixed)
        layout.addWidget(line, 4, 0, 2, 4)

        # input parameters

        self.export_input_parameters_QPushButton = QPushButton("Export input (geographical-geological) parameters")
        self.export_input_parameters_QPushButton.clicked.connect(self.export_input_parameters)
        layout.addWidget(self.export_input_parameters_QPushButton, 6, 0, 1, 1)

        self.export_params_file_path_QLineEdit = QLineEdit()
        layout.addWidget(self.export_params_file_path_QLineEdit, 6, 1, 1, 2)

        self.define_export_file_path_QPushButton = QPushButton(".....")
        self.define_export_file_path_QPushButton.clicked.connect(self.define_output_params_file_path)
        layout.addWidget(self.define_export_file_path_QPushButton, 6, 3, 1, 1)

        widget.setLayout(layout)

        return widget

    def export_input_parameters(self):

        if not hasattr(self, 'source_point_longitude'):
            QMessageBox.critical(
                self,
                self.module_name,
                f"Source point longitude is undefined"
            )
            return

        if not hasattr(self, 'source_point_latitude'):
            QMessageBox.critical(
                self,
                self.module_name,
                f"Source point latitude is undefined"
            )
            return

        if not hasattr(self, 'source_point_elevation'):
            QMessageBox.critical(
                self,
                self.module_name,
                f"Source point elevation is undefined"
            )
            return

        if not hasattr(self, 'dip_direction_as_north_azimuth'):
            QMessageBox.critical(
                self,
                self.module_name,
                f"Dip direction (relative to Geographic North) not calculated"
            )
            return

        if not hasattr(self, 'source_dip_angle'):
            QMessageBox.critical(
                self,
                self.module_name,
                f"Dip angle is undefined"
            )
            return

        export_file_path = self.export_params_file_path_QLineEdit.text()

        if export_file_path == "":
            QMessageBox.critical(
                self,
                self.module_name,
                "No export file path defined"
            )
            return

        with open(export_file_path, "w") as ofile:
            ofile.write(f"longitude: {self.source_point_longitude}\n")
            ofile.write(f"latitude: {self.source_point_latitude}\n")
            ofile.write(f"elevation: {self.source_point_elevation}\n")
            ofile.write(f"dip direction azimuth (geographic North): {self.dip_direction_as_north_azimuth}\n")
            ofile.write(f"dip angle: {self.source_dip_angle}")

        QMessageBox.information(
            self,
            self.module_name,
            f"Processing parameters exported as file '{export_file_path}'"
        )
        return

    def load_input_processing_parameters(self):

        try:

            input_file_path, _ = QFileDialog.getOpenFileName(
                self,
                self.tr("Load processing parameters"),
                "*.txt",
                "txt (*.txt *.TXT)"
            )

            if not input_file_path:
                QMessageBox.critical(
                    self,
                    self.module_name,
                    "No input file path is defined"
                )
                return

            with open(input_file_path, "r") as ifile:
                raw_params = ifile.readlines()

            params = list(map(lambda val: val.split(":")[-1].strip(), raw_params))

            self.lock_z_value_to_dem_surface_QCheckBox.setChecked(False)

            self.source_point_longitude_QLineEdit.setText(params[0])
            self.source_point_longitude = float(params[0])
            self.source_point_latitude_QLineEdit.setText(params[1])
            self.source_point_latitude = float(params[1])
            self.source_point_elevation_QLineEdit.setText(params[2])
            self.source_point_elevation = float(params[2])

            self.direction_azimuth_is_relative_to_geographic_north_QRadioButton.setChecked(True)
            self.dip_direction_QDoubleSpinBox.setValue(float(params[3]))
            self.dip_direction_as_north_azimuth = float(params[3])
            self.dip_angle_QDoubleSpinBox.setValue(float(params[4]))
            self.source_dip_angle = float(params[4])

            self.update_source_point_location()

            QMessageBox.information(
                self,
                self.module_name,
                f"Parameters loaded"
            )
            return

        except Exception as e:

            QMessageBox.critical(
                self,
                self.module_name,
                f"Error while reading input file:\n{e!r}"
            )
            return

    def set_source_point_in_map(self):
        
        try:
            self.map_tool.canvasClicked.disconnect(self.update_source_point_location)
        except:
            pass            
         
        self.update_crs_settings()

        # mouse listener
        self.map_tool = EmitPointFromCanvas(
            self.canvas,
            self.plugin
        )

        # save the standard map tool for restoring it at the end
        self.previous_tool = self.canvas.mapTool()
        self.map_tool.canvasClicked.connect(self.update_source_point_location)
        self.map_tool.setCursor(Qt.CursorShape.CrossCursor)
        self.canvas.setMapTool(self.map_tool)

    def get_chosen_dem(self, ndx_dem_file = 0):

        self.dem = None        
        self.reset_post_dem_change()
        if self.raster_layers is None or len(self.raster_layers) == 0:
            return          
                                
        # no DEM layer defined  
        if ndx_dem_file == 0:
            return             

        self.dem = self.raster_layers[ndx_dem_file - 1]

        dem_src = self.dem.source()

        results, err = read_raster_band_with_wkt(
            raster_source=dem_src)

        if err:
            QMessageBox.critical(
                self,
                self.module_name,
                f"Exception while reading DEM file:\n{err!r}"
            )
            return

        geotransform, band_params, data, wkt = results

        self.grid = Grid(
            array=data,
            geotransform=geotransform,
            wkt_crs=wkt,
        )

        self.set_source_point_in_map_QPushButton.setEnabled(True)
        self.reset_source_point_in_map_QPushButton.setEnabled(True)

    def create_marker(
        self,
            canvas,
            prj_crs_x,
            prj_crs_y,
            pen_width=UIConstants.DEFAULT_PEN_WIDTH,
            icon_type=UIConstants.DEFAULT_ICON_TYPE,
            icon_size=UIConstants.DEFAULT_ICON_SIZE,
            icon_color: Union[str, QColor] = UIConstants.DEFAULT_ICON_COLOR
    ) -> QgsVertexMarker:
        
        marker = QgsVertexMarker(canvas)
        marker.setIconType(icon_type)
        marker.setIconSize(icon_size)
        marker.setPenWidth(pen_width)
        if isinstance(icon_color, QColor):
            color = icon_color
        else:
            color = QColor(icon_color)
        marker.setColor(color)
        marker.setCenter(QgsPointXY(prj_crs_x, prj_crs_y))

        return marker        

    def update_marker_in_canvas(self) -> bool:

        try:

            self.remove_intersection_markers_from_canvas()
            self.remove_source_point_marker_from_canvas()

            if self.source_point_longitude is not None and isfinite(self.source_point_longitude) and \
                    self.source_point_latitude is not None and isfinite(self.source_point_latitude):
                color = QColor(str(self.intersections_color_QComboBox.currentText()))

                result, err = qgis_project_coords(
                    x=self.source_point_longitude,
                    y=self.source_point_latitude,
                    dest_crs=self.map_crs,
                )

                if err:

                    QMessageBox.critical(
                        self,
                        self.module_name,
                        f"{err!r}"
                    )

                    return True

                source_pt_x, source_pt_y = result

                source_point_marker = self.create_marker(
                    self.canvas,
                    source_pt_x,
                    source_pt_y,
                    icon_type=1,
                    icon_color=color)

                self.source_point_marker_list = [source_point_marker]

            self.canvas.refresh()

            return False

        except Exception as e:

            QMessageBox.critical(
                self,
                self.module_name,
                f"Exception while updating marker:\n{e!r}"
            )

            return True

    def update_source_point_location(
        self,
        qgs_point=None,
    ):

        try:

            self.reset_source_point_in_map_QPushButton.setEnabled(True)

            if self.sender() == self.map_tool:

                result, err = qgis_project_coords(
                    x=qgs_point.x(),
                    y=qgs_point.y(),
                    src_crs=self.map_crs,
                )

                if err:
                    QMessageBox.critical(
                        self,
                        self.module_name,
                        f"{err!r}"
                    )
                    return

                self.source_point_longitude, self.source_point_latitude = result

                self.source_point_longitude_QLineEdit.setText(str(self.source_point_longitude))
                self.source_point_latitude_QLineEdit.setText(str(self.source_point_latitude))

            elif self.sender() == self.source_point_longitude_QLineEdit:

                if self.source_point_longitude_QLineEdit.text() == '':
                    QMessageBox.critical(
                        self,
                        self.module_name,
                        "Error: longitude value is not defined"
                    )
                    return

                try:

                    self.source_point_longitude = float(self.source_point_longitude_QLineEdit.text())

                except Exception as e:

                    QMessageBox.critical(
                        self,
                        self.module_name,
                        f"Longitude value is not correctly defined:\n{e!r}"
                    )
                    return

            elif self.sender() == self.source_point_latitude_QLineEdit:

                if self.source_point_latitude_QLineEdit.text() == '':
                    QMessageBox.critical(
                        self,
                        self.module_name,
                        "Latitude value is undefined"
                    )
                    return

                try:

                    self.source_point_latitude = float(self.source_point_latitude_QLineEdit.text())

                except Exception as e:
                    QMessageBox.critical(
                        self,
                        self.module_name,
                        f"Latitude value is not correctly defined:\n{e!r}"
                    )
                    return

            if self.lock_z_value_to_dem_surface_QCheckBox.isChecked():

                elevation_value_from_dem = self.update_source_point_elevation_value()

                if elevation_value_from_dem is None:
                    self.current_source_point_elevation_value = None
                    self.source_point_elevation_QLineEdit.setText("")
                else:
                    self.current_source_point_elevation_value = elevation_value_from_dem
                    self.source_point_elevation_QLineEdit.setText(str(self.current_source_point_elevation_value))
                    self.use_elevation_from_dem = True

            _ = self.update_marker_in_canvas()

        except Exception as e:

            QMessageBox.critical(
                self,
                self.module_name,
                f"Exception while updating source point location:\n{e!r}"
            )
            return

    def update_source_point_elevation_value (self):
        """
        Update z value.
        
        """

        try:

            # does nothing when the height source is not from DEM

            if not self.lock_z_value_to_dem_surface_QCheckBox.isChecked():
                return None

            # prevent action when the DEM is not read

            if self.grid is None:
                return None

            # get z value from ptLyr

            if self.source_point_longitude is None or self.source_point_latitude is None:
                return None

            result, err = qgis_project_coords(
                x=self.source_point_longitude,
                y=self.source_point_latitude,
                dest_crs=self.dem.crs(),
            )

            if err:
                QMessageBox.critical(
                    self,
                    self.module_name,
                    f"{err!r}"
                )
                return

            srcpt_demcrs_x, srcpt_demcrs_y = result

            # return None or a numeric value
            return self.grid.interpolate_bilinear(
                srcpt_demcrs_x,
                srcpt_demcrs_y)

        except Exception as e:

            QMessageBox.critical(
                self,
                self.module_name,
                f"Exception while updating z value:\n{e!r}"
            )

            return

    def restore_previous_map_tool(self,
                                  map_tool):

        if map_tool:
            self.canvas.unsetMapTool(map_tool)

        if self.previous_tool:
            self.canvas.setMapTool(self.previous_tool)

    def clear_source_point_qlineedit_forms(self):
        
        for qlineedit_form in (
                self.source_point_longitude_QLineEdit,
                self.source_point_latitude_QLineEdit,
                self.source_point_elevation_QLineEdit):

            qlineedit_form.clear()

    def update_dip_direction_slider(self):
        """
        Update the value of the dip direction in the slider.""
        """        
        transformed_dipdirection = self.dip_direction_QDoubleSpinBox.value() + 180.0
        if transformed_dipdirection > GeometryConstants.DIP_DIRECTION_NORMALIZED_MAX:
            transformed_dipdirection -= GeometryConstants.DIP_DIRECTION_NORMALIZED_MAX
        self.dip_direction_QDial.setValue(int(transformed_dipdirection))

    def update_dip_direction_spinbox(self):
        """
        Update the value of the dip direction in the spinbox.
        """        
        real_dipdirection = self.dip_direction_QDial.value() - 180.0
        if real_dipdirection < GeometryConstants.DIP_DIRECTION_NORMALIZED_MIN:
            real_dipdirection += GeometryConstants.DIP_DIRECTION_NORMALIZED_MAX
        self.dip_direction_QDoubleSpinBox.setValue(real_dipdirection)

    def update_dip_angle_slider(self):
        """
        Update the value of the dip angle in the slider.
        """
        
        self.DAngle_verticalSlider.setValue(int(self.dip_angle_QDoubleSpinBox.value()))

    def update_dip_angle_spinbox(self):
        """
        Update the value of the dip angle in the spinbox.
        """        
        
        self.dip_angle_QDoubleSpinBox.setValue(self.DAngle_verticalSlider.value())

    def project_from_dem_to_map_crs(
        self,
        x: float,
        y: float
    ) -> Tuple[Union[None, Tuple[float, float], Error]]:
        """
        Converts point coordinates from project to DEM crs.

        :param x: x coordinate in the DEM CRS.
        :param y: y coordinate in the DEM CRS.
        :return: the point coordinates in the project crs.
        """

        return qgis_project_coords(
            x,
            y,
            self.dem.crs(),
            self.map_crs
        )

    def parse_input_parameters(self) -> Tuple[Union[None, numbers.Real], Error]:

        try:

            # check if all input data are correct

            if self.grid is None:
                z = None
                msg = "Please first define a source DEM"
            elif self.source_point_longitude_QLineEdit.text() == '' or self.source_point_latitude_QLineEdit.text() == '' or self.source_point_elevation_QLineEdit.text() == '':
                z = None
                msg = "Define the location of the source point in 'Geographic parameters' section"
            elif self.source_point_elevation_QLineEdit.text() == '':
                z = None
                msg = "Define the elevation of the source point in 'Geographic parameters' section"
            else:
                try:
                    z = float(self.source_point_elevation_QLineEdit.text())
                    msg = ''
                except:
                    z = None
                    msg = "z value is not correctly defined"

            if z is None:

                err = Error(
                    True,
                    caller_name(),
                    Exception(msg),
                    traceback.format_exc(),
                )

            else:

                err = Error()

            return z, err

        except Exception as e:

            return None, Error(
                True,
                caller_name(),
                e,
                traceback.format_exc(),
            )

    def preprocess_plane_and_source_point(
        self,
        local_pt_longitude: numbers.Real,
        local_pt_latitude: numbers.Real,
    ) -> Tuple[Union[None, Tuple[Plane, Point]], Error]:

        try:

            self.source_dip_direction = self.dip_direction_QDoubleSpinBox.value()
            self.source_dip_angle = self.dip_angle_QDoubleSpinBox.value()

            # Calculates dip direction correction with respect to map CRS top (y-axis) direction

            if self.direction_azimuth_is_relative_to_map_top_QRadioButton.isChecked():

                geographic_north_direction_azimuth_cw, err = calculate_local_geographic_north_azimuth_from_epsg_4326(
                    longitude=local_pt_longitude,
                    latitude=local_pt_latitude,
                    planar_qgis_crs=self.map_crs,
                ) #-> Tuple[Union[None, numbers.Real], Error]

                if err:
                    return None, err

                self.dip_direction_as_north_azimuth = (self.source_dip_direction - geographic_north_direction_azimuth_cw) % 360.0

            else:

                self.dip_direction_as_north_azimuth = self.source_dip_direction % 360.0

            # Correct dip direction and angle by DEM crs

            dem_crs_plane, err = convert_plane_orientation_to_crs(
                dip_direction_as_north_azimuth=self.dip_direction_as_north_azimuth,
                source_dip_angle=self.source_dip_angle,
                local_pt_longitude=local_pt_longitude,
                local_pt_latitude=local_pt_latitude,
                dest_crs=self.dem.crs(),
            )

            if err:
                return None, err

            result, err = qgis_project_coords(
                x=local_pt_longitude,
                y=local_pt_latitude,
                dest_crs=self.dem.crs(),
            )

            if err:
                return None, err

            dem_crs_src_pt_x, dem_crs_src_pt_y = result

            dem_crs_src_pt_3d = Point(
                dem_crs_src_pt_x,
                dem_crs_src_pt_y,
                self.source_point_elevation,
            )

            return (
                        dem_crs_plane,
                        dem_crs_src_pt_3d
                    ), Error()

        except Exception as e:

            return None, Error(
                True,
                caller_name(),
                e,
                traceback.format_exc(),
            )

    def calculate_intersections(self):
        """
        Calcola punti di intersezione tra piano geologico e DEM.
        Coordina parsing parametri, preprocessing, calcolo intersezioni e preparazione risultati.
        """
        try:
            # Parsing parametri input
            elevation, err = self.parse_input_parameters()
            if err:
                QMessageBox.critical(
                    self,
                    self.module_name,
                    f"Exception while parsing input parameters for intersection calculation:\n{err!r}"
                )
                return

            self.source_point_elevation = elevation

            # Preprocessing piano e punto sorgente
            plane_and_point, err = self._preprocess_plane_and_point()
            if err:
                QMessageBox.critical(
                    self,
                    self.module_name,
                    f"Exception while preprocessing plane and source point:\n{err!r}"
                )
                return

            self.geological_plane_in_dem_crs, source_point_in_dem_crs = plane_and_point

            # Calcolo intersezioni
            intersection_coords, err = self._calculate_plane_dem_intersections(
                source_point_in_dem_crs
            )
            if err:
                QMessageBox.critical(
                    self,
                    self.module_name,
                    f"Exception while calculating DEM-plane intersections:\n{err!r}"
                )
                return

            # Preparazione risultati
            self._prepare_intersection_results(intersection_coords, source_point_in_dem_crs)
            self.plot_intersections()

        except Exception as e:
            QMessageBox.critical(
                self,
                self.module_name,
                f"Exception while calculating intersections:\n{e!r}"
            )

    def _preprocess_plane_and_point(self):
        """Preprocessa piano e punto sorgente."""
        return self.preprocess_plane_and_source_point(
            self.source_point_longitude,
            self.source_point_latitude
        )

    def _calculate_plane_dem_intersections(self, source_point_in_dem_crs):
        """Calcola intersezioni piano-DEM con o senza buffer."""
        if self.use_buffer_zone_QCheckBox.isChecked():
            buffer_radius = self.buffer_radius_QDoubleSpinBox.value()
            result = plane_dem_intersections_with_buffer(
                plane_attitude=self.geological_plane_in_dem_crs,
                source_point=source_point_in_dem_crs,
                grid=self.grid,
                buffer_radius=buffer_radius
            )
        else:
            result = plane_dem_intersections_strategy_1(
                plane_attitude=self.geological_plane_in_dem_crs,
                source_point=source_point_in_dem_crs,
                grid=self.grid
            )

        if isinstance(result, Error):
            return None, result

        return result, None

    def _prepare_intersection_results(self, intersection_coords, source_point_in_dem_crs):
        """Prepara risultati intersezioni per visualizzazione ed export."""
        xcoords_x, xcoords_y, ycoords_x, ycoords_y = intersection_coords

        # Crea oggetto Intersections
        self.inters = Intersections()
        self.inters.xcoords_x = xcoords_x
        self.inters.xcoords_y = xcoords_y
        self.inters.ycoords_x = ycoords_x
        self.inters.ycoords_y = ycoords_y
        self.inters.parameters = Intersection_Parameters(
            source_point_in_dem_crs,
            self.geological_plane_in_dem_crs
        )

        # Estrae coordinate valide (non NaN)
        self.intersections_x = list(xcoords_x[np.logical_not(np.isnan(xcoords_x))]) + \
                               list(ycoords_x[np.logical_not(np.isnan(ycoords_y))])
        self.intersections_y = list(xcoords_y[np.logical_not(np.isnan(xcoords_x))]) + \
                               list(ycoords_y[np.logical_not(np.isnan(ycoords_y))])

        # Prepara dati per export
        intersection_data = dict(x=self.intersections_x, y=self.intersections_y)
        intersection_plane = dict(
            dipdir=self.inters.parameters._srcPlaneAttitude._dipdir,
            dipangle=self.inters.parameters._srcPlaneAttitude._dipangle
        )
        intersection_point = dict(
            x=self.inters.parameters._srcPt.x,
            y=self.inters.parameters._srcPt.y,
            z=self.inters.parameters._srcPt.z
        )

        self.intersections_xprt = dict(
            data=intersection_data,
            plane=intersection_plane,
            point=intersection_point
        )

    def plot_intersections(self):

        try:

            current_markers_list = []

            for x, y in zip(self.intersections_x, self.intersections_y):

                result, err = self.project_from_dem_to_map_crs(x, y)
                if err:
                    return err

                prj_crs_x, prj_crs_y = result

                marker = self.create_marker(
                    self.canvas,
                    prj_crs_x,
                    prj_crs_y,
                    pen_width=1,
                    icon_type=1,
                    icon_size=8,
                    icon_color=QColor(str(self.intersections_color_QComboBox.currentText())))

                current_markers_list.append(marker)

            self.intersections_markers_list += current_markers_list
            self.canvas.refresh()

        except Exception as e:

            return Error(
                True,
                caller_name(),
                e,
                traceback.format_exc(),
            )

    def define_output_params_file_path(self):

        output_filename, _ = QFileDialog.getSaveFileName(
            self,
            self.tr("Save processing parameters"),
            "*.txt",
            "txt (*.txt *.TXT)"
        )

        if not output_filename:
            return

        self.export_params_file_path_QLineEdit.setText(output_filename)

    def define_output_vector_file_path(self):
            
        output_filename, _ = QFileDialog.getSaveFileName(
            self,
            self.tr("Save shapefile"),
            "*.shp",
            "shp (*.shp *.SHP)"
        )

        if not output_filename: 
            return

        self.output_shapefile_path_QLineEdit.setText(output_filename)

    def load_output_grid_parameters(self):

        output_filename, __ = QFileDialog.getSaveFileName(
            self,
            self.tr("Save grid as geotiff"),
            "*.tif",
            "tif (*.tif *.TIF)"
        )

        if not output_filename:
            return

        self.output_grid_path_input_QLineEdit.setText(output_filename)

    def select_output_grid_file(self):

        output_filename, __ = QFileDialog.getSaveFileName(
            self,
            self.tr("Save grid as geotiff"),
            "*.tif",
            "tif (*.tif *.TIF)"
        )

        if not output_filename:
            return

        self.output_grid_path_input_QLineEdit.setText(output_filename)

    def export_last_intersection_result(self):
        """
        Write intersection results in the output shapefile.
        """

        # check for result existence

        if self.inters is None:
            QMessageBox.critical(
                self,
                self.module_name,
                "No results available"
            )
            return

        self.output_filename = str(self.output_shapefile_path_QLineEdit.text())

        if self.output_filename == '':
            QMessageBox.critical(
                self,
                self.module_name,
                "No output file is defined"
            )
            return

        # set output type

        if self.save_intersections_as_points_QRadioButton.isChecked():

            self.result_geometry = 'points'

        else:

            self.result_geometry = 'lines'

        # creation of output shapefile

        shape_driver = ogr.GetDriverByName("ESRI Shapefile")

        self.out_shape = shape_driver.CreateDataSource(
            self.output_filename,
        )

        if self.out_shape is None:
            QMessageBox.critical(
                self,
                self.module_name,
                f"Unable to create output shapefile '{self.output_filename}'"
            )
            return

        dest_srs = ogr.osr.SpatialReference()
        dest_srs.ImportFromEPSG(EPSGCodes.WGS84)
        # Use 3D geometry types since we always have z values
        if self.result_geometry == 'points':
            geom_type = ogr.wkbPoint25D
        else:
            geom_type = ogr.wkbLineString25D

        self.out_layer = self.out_shape.CreateLayer(
            'output_lines',
            dest_srs,
            geom_type=geom_type,
        )

        if self.out_layer is None:
            QMessageBox.critical(
                self,
                self.module_name,
                f"Unable to create output shapefile '{self.output_filename}'",
            )
            return

        # set analysis parameters
        sourcePoint = self.inters.parameters._srcPt
        self.srcPlaneAttitude = self.inters.parameters._srcPlaneAttitude
        self.z_from_plane_closure = self.srcPlaneAttitude.plane_from_geo(sourcePoint)

        # add fields to the output shapefile
        id_fieldDef = ogr.FieldDefn(ID_FLD_NM, ogr.OFTInteger)
        self.out_layer.CreateField(id_fieldDef)

        if self.result_geometry == 'points':
            longitude_fieldDef = ogr.FieldDefn(LON_FLD_NM, ogr.OFTReal)
            self.out_layer.CreateField(longitude_fieldDef)
            latitude_fieldDef = ogr.FieldDefn(LAT_FLD_NM, ogr.OFTReal)
            self.out_layer.CreateField(latitude_fieldDef)
            z_fieldDef = ogr.FieldDefn(Z_FLD_NM, ogr.OFTReal)
            self.out_layer.CreateField(z_fieldDef)

        src_pt_lon_field_def = ogr.FieldDefn(SRC_PT_LON_FLD_NM, ogr.OFTReal)
        self.out_layer.CreateField(src_pt_lon_field_def)
        src_pt_lat_field_def = ogr.FieldDefn(SRC_PT_LAT_FLD_NM, ogr.OFTReal)
        self.out_layer.CreateField(src_pt_lat_field_def)
        src_pt_z_field_def = ogr.FieldDefn(SRC_PT_Z_FLD_NM, ogr.OFTReal)
        self.out_layer.CreateField(src_pt_z_field_def)
        dip_dir_field_def = ogr.FieldDefn(DIP_DIR_FLD_NM, ogr.OFTReal)
        self.out_layer.CreateField(dip_dir_field_def)
        dip_ang_field_def = ogr.FieldDefn(DIP_ANG_FLD_NM, ogr.OFTReal)
        self.out_layer.CreateField(dip_ang_field_def)

        # get the layer definition of the output shapefile
        self.output_shape_feature_definitions = self.out_layer.GetLayerDefn()

        # write results
        if self.result_geometry == 'points':
            self.write_intersections_as_points()
        else:
            self.write_intersections_as_lines()

        QMessageBox.information(
            self,
            self.module_name,
            f"Saved as shapefile '{self.output_filename}'"
        )

        # add theme to QGis project
        if self.load_output_in_map_QCheckBox.isChecked():
            try:
                intersection_layer = QgsVectorLayer(
                    self.output_filename,
                    QFileInfo(self.output_filename).baseName(),
                    "ogr")

                # Verify layer is valid before adding to project
                if not intersection_layer.isValid():
                    QMessageBox.critical(
                        self,
                        self.module_name,
                        f"Unable to load layer '{intersection_layer.name()}' from '{self.output_filename}': layer is not valid"
                    )
                    return

                QgsProject.instance().addMapLayer(intersection_layer)
            except Exception as e:
                QMessageBox.critical(
                    self,
                    self.module_name,
                    f"Unable to load layer in project:\n{e!r}"
                )
                return

    def preprocess_intersections_for_export_as_points(self
    ) -> Tuple[Union[None, Tuple[numbers.Real, numbers.Real, numbers.Real,]], Error]:

        try:

            x_filtered_coord_x = self.inters.xcoords_x[np.logical_not(np.isnan(self.inters.xcoords_x))]
            x_filtered_coord_y = self.inters.xcoords_y[np.logical_not(np.isnan(self.inters.xcoords_x))]
            x_filtered_coord_z = self.z_from_plane_closure(x_filtered_coord_x, x_filtered_coord_y)

            y_filtered_coord_x = self.inters.ycoords_x[np.logical_not(np.isnan(self.inters.ycoords_y))]
            y_filtered_coord_y = self.inters.ycoords_y[np.logical_not(np.isnan(self.inters.ycoords_y))]
            y_filtered_coord_z = self.z_from_plane_closure(y_filtered_coord_x, y_filtered_coord_y)

            intersections_x = list(x_filtered_coord_x) + list(y_filtered_coord_x)
            intersections_y = list(x_filtered_coord_y) + list(y_filtered_coord_y)
            intersections_z = list(x_filtered_coord_z) + list(y_filtered_coord_z)

            points_in_epsg_4326 = []

            for dem_crs_x, dem_crs_y, z in zip(intersections_x, intersections_y, intersections_z):

                point_coords, err = qgis_project_coords(
                    x=dem_crs_x,
                    y=dem_crs_y,
                    src_crs=self.dem.crs(),
                )

                if err:
                    QMessageBox.critical(
                        self,
                        self.module_name,
                        f"Error while preprocessing intersection coordinates for export in point shapefile:\n{err!r}"
                    )
                    return

                pt_longitude, pt_latitude = point_coords

                points_in_epsg_4326.append((pt_longitude, pt_latitude, z))

            return points_in_epsg_4326, Error()

        except Exception as e:

            return None, Error(
                True,
                caller_name(),
                e,
                traceback.format_exc(),
            )

    def write_intersections_as_points(self):
        """
        Write intersection results in the output point shapefile.
        """

        points_in_epsg_4326, err = self.preprocess_intersections_for_export_as_points()

        if err:
            QMessageBox.critical(
                self,
                self.module_name,
                f"Error while preprocessing intersection coordinates for export in point shapefile:\n{err!r}"
            )
            return

        curr_pt_id = 0

        for lon, lat, z in points_in_epsg_4326:

            curr_pt_id += 1

            # Validate and convert coordinates to float
            try:
                lon_f = float(lon)
                lat_f = float(lat)
                z_f = float(z) if z is not None else None
            except (TypeError, ValueError) as e:
                QgsMessageLog.logMessage(
                    f"Invalid coordinate values: lon={lon}, lat={lat}, z={z}, error: {e!r}",
                    "qgSurf",
                    Qgis.MessageLevel.Warning
                )
                continue

            # Validate that coordinates are finite
            if not (math.isfinite(lon_f) and math.isfinite(lat_f)):
                QgsMessageLog.logMessage(
                    f"Non-finite coordinates: lon={lon_f}, lat={lat_f}",
                    "qgSurf",
                    Qgis.MessageLevel.Warning
                )
                continue

            # Determine geometry type based on z value
            if z_f is not None and math.isfinite(z_f):
                curr_pt_geom = ogr.Geometry(ogr.wkbPoint25D)
                curr_pt_geom.AddPoint(lon_f, lat_f, z_f)
            else:
                curr_pt_geom = ogr.Geometry(ogr.wkbPoint)
                curr_pt_geom.AddPoint(lon_f, lat_f)

            # create a new feature
            curr_Pt_shape = ogr.Feature(self.output_shape_feature_definitions)
            curr_Pt_shape.SetGeometry(curr_pt_geom)
            curr_Pt_shape.SetField(ID_FLD_NM, curr_pt_id)

            curr_Pt_shape.SetField(LON_FLD_NM, lon_f)
            curr_Pt_shape.SetField(LAT_FLD_NM, lat_f)
            if z_f is not None and math.isfinite(z_f):
                curr_Pt_shape.SetField(Z_FLD_NM, z_f)

            curr_Pt_shape.SetField(SRC_PT_LON_FLD_NM, self.source_point_longitude)
            curr_Pt_shape.SetField(SRC_PT_LAT_FLD_NM, self.source_point_latitude)
            curr_Pt_shape.SetField(SRC_PT_Z_FLD_NM, self.source_point_elevation)

            curr_Pt_shape.SetField(DIP_DIR_FLD_NM, self.dip_direction_as_north_azimuth)
            curr_Pt_shape.SetField(DIP_ANG_FLD_NM, self.source_dip_angle)

            # add the feature to the output layer and verify success
            result = self.out_layer.CreateFeature(curr_Pt_shape)
            if result != ogr.OGRERR_NONE:
                QgsMessageLog.logMessage(
                    f"Error creating feature {curr_pt_id}: OGR error code {result}",
                    "qgSurf",
                    Qgis.MessageLevel.Warning
                )

            # destroy no longer used objects
            try:
                curr_pt_geom.Destroy()
            except:
                pass
            try:
                curr_Pt_shape.Destroy()
            except:
                pass

        # destroy output geometry
        try:
            self.out_shape.Destroy()
        except:
            pass

    def write_intersections_as_lines(self):
        """
        Write intersection results in a line shapefile.
        """

        # create dictionary of cell with intersection points
        self.inters.links = self.inters.get_intersections()
        self.inters.neighbours = self.inters.set_neighbours()
        self.inters.define_paths()

        # networks of connected intersections
        self.inters.networks = self.inters.define_networks()

        for curr_path_id, curr_path_points in self.inters.networks.items():

            # Determine if we need 3D geometry by checking first point
            first_point_checked = False
            needs_3d = False
            
            # First pass: determine if we need 3D geometry
            for curr_point_id in curr_path_points:
                curr_intersection = self.inters.links[curr_point_id - 1]
                i, j, direct = curr_intersection['i'], curr_intersection['j'], curr_intersection['pi_dir']
                if direct == 'x':
                    dem_crs_x, dem_crs_y = self.inters.xcoords_x[i, j], self.inters.xcoords_y[i, j]
                elif direct == 'y':
                    dem_crs_x, dem_crs_y = self.inters.ycoords_x[i, j], self.inters.ycoords_y[i, j]
                else:
                    continue
                    
                z = self.z_from_plane_closure(dem_crs_x, dem_crs_y)
                if z is not None and math.isfinite(z):
                    needs_3d = True
                break  # Check only first point
            
            # Create geometry with appropriate type
            if needs_3d:
                line = ogr.Geometry(ogr.wkbLineString25D)
            else:
                line = ogr.Geometry(ogr.wkbLineString)

            # Second pass: add points
            for curr_point_id in curr_path_points:
                curr_intersection = self.inters.links[curr_point_id - 1]
                i, j, direct = curr_intersection['i'], curr_intersection['j'], curr_intersection['pi_dir']
                if direct == 'x':
                    dem_crs_x, dem_crs_y = self.inters.xcoords_x[i, j], self.inters.xcoords_y[i, j]
                elif direct == 'y':
                    dem_crs_x, dem_crs_y = self.inters.ycoords_x[i, j], self.inters.ycoords_y[i, j]
                else:
                    continue
                    
                z = self.z_from_plane_closure(dem_crs_x, dem_crs_y)

                point_coords, err = qgis_project_coords(
                    x=dem_crs_x,
                    y=dem_crs_y,
                    src_crs=self.dem.crs(),
                )

                if err:
                    QMessageBox.critical(
                        self,
                        self.module_name,
                        f"Error while exporting coordinates in line shapefile:\n{err!r}"
                    )
                    return

                pt_longitude, pt_latitude = point_coords
                
                # Validate and convert to float
                try:
                    lon_f = float(pt_longitude)
                    lat_f = float(pt_latitude)
                    z_f = float(z) if z is not None and math.isfinite(z) else None
                except (TypeError, ValueError) as e:
                    QgsMessageLog.logMessage(
                        f"Invalid coordinate values: lon={pt_longitude}, lat={pt_latitude}, z={z}, error: {e!r}",
                        "qgSurf",
                        Qgis.MessageLevel.Warning
                    )
                    continue

                # Validate that coordinates are finite
                if not (math.isfinite(lon_f) and math.isfinite(lat_f)):
                    QgsMessageLog.logMessage(
                        f"Non-finite coordinates: lon={lon_f}, lat={lat_f}",
                        "qgSurf",
                        Qgis.MessageLevel.Warning
                    )
                    continue

                # Add point with appropriate number of coordinates
                if needs_3d and z_f is not None and math.isfinite(z_f):
                    line.AddPoint(lon_f, lat_f, z_f)
                else:
                    line.AddPoint(lon_f, lat_f)

            # Skip empty geometries
            if line.GetPointCount() == 0:
                QgsMessageLog.logMessage(
                    f"Skipping empty line geometry for path {curr_path_id}",
                    "qgSurf",
                    Qgis.MessageLevel.Warning
                )
                continue

            # create a new feature
            line_shape = ogr.Feature(self.output_shape_feature_definitions)
            line_shape.SetGeometry(line)

            line_shape.SetField(ID_FLD_NM, curr_path_id)

            line_shape.SetField(SRC_PT_LON_FLD_NM, self.source_point_longitude)
            line_shape.SetField(SRC_PT_LAT_FLD_NM, self.source_point_latitude)
            line_shape.SetField(SRC_PT_Z_FLD_NM, self.source_point_elevation)

            line_shape.SetField(DIP_DIR_FLD_NM, self.dip_direction_as_north_azimuth)
            line_shape.SetField(DIP_ANG_FLD_NM, self.source_dip_angle)

            # add the feature to the output layer and verify success
            result = self.out_layer.CreateFeature(line_shape)
            if result != ogr.OGRERR_NONE:
                QgsMessageLog.logMessage(
                    f"Error creating line feature {curr_path_id}: OGR error code {result}",
                    "qgSurf",
                    Qgis.MessageLevel.Warning
                )

            # destroy no longer used objects
            try:
                line_shape.Destroy()
            except:
                pass

        # destroy output geometry
        self.out_shape.Destroy()


