# -*- coding: utf-8 -*-
"""
Dual Profile Tool - Main dialog and drawing tool
Adapted from the original script for plugin structure
"""

from qgis.PyQt import QtWidgets, QtCore
from qgis.PyQt.QtGui import QColor, QFont
from qgis.PyQt.QtCore import Qt, QVariant
from . import qt_compat

from qgis.core import (
    # Basics
    QgsProject, QgsVectorLayer, QgsFeature, QgsGeometry,
    QgsPoint, QgsPointXY, QgsLineString, QgsWkbTypes,
    QgsField, QgsFields, QgsVectorFileWriter,
    QgsCoordinateTransform, QgsCoordinateReferenceSystem,
    
    # Raster
    QgsRasterLayer, QgsRaster,
    
    # Symbology
    QgsSymbol, QgsSimpleLineSymbolLayer,
    QgsSingleSymbolRenderer, QgsMarkerSymbol,
    QgsSimpleMarkerSymbolLayer,
    QgsCategorizedSymbolRenderer, QgsRendererCategory,
    
    # Labeling
    QgsTextFormat, QgsVectorLayerSimpleLabeling,
    QgsPalLayerSettings, QgsTextBufferSettings,
    
    # General
    Qgis
)
from qgis.gui import QgsMapTool, QgsRubberBand, QgsMapToolEmitPoint
import math
import numpy as np
from datetime import datetime
import os
import tempfile
import webbrowser

# Import profile exporter
from .profile_exporter import ProfileExporter
from .vector_export_dialog import VectorExportDialog
from .multi_dem_widget import MultiDEMWidget

# Try to import plotly
try:
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots
    import plotly.offline as pyo
    PLOTLY_AVAILABLE = True
except ImportError:
    PLOTLY_AVAILABLE = False
    print("⚠️ plotly not available. Install with: pip install plotly")

# Fallback to matplotlib
if not PLOTLY_AVAILABLE:
    try:
        import matplotlib
        matplotlib.use('Qt5Agg')
        from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
        from matplotlib.figure import Figure
        import matplotlib.pyplot as plt
        MATPLOTLIB_AVAILABLE = True
    except ImportError:
        MATPLOTLIB_AVAILABLE = False

class DualProfileTool(QgsMapTool):
    """Tool for drawing two parallel lines centered on the drawn line"""
    
    def __init__(self, canvas, callback):
        super().__init__(canvas)
        self.canvas = canvas
        self.callback = callback
        self.rubber_band_center = None
        self.rubber_band1 = None
        self.rubber_band2 = None
        self.start_point = None
        self.end_point = None
        self.is_drawing = False
        self.offset_distance = 1.0
        
    def canvasPressEvent(self, event):
        """Handle mouse clicks"""
        if event.button() == Qt.LeftButton:
            point = self.toMapCoordinates(event.pos())
            
            if not self.is_drawing:
                # First click: start drawing
                self.start_point = point
                self.is_drawing = True
                
                # Create rubber bands
                self.rubber_band_center = QgsRubberBand(self.canvas, QgsWkbTypes.LineGeometry)
                self.rubber_band_center.setColor(QColor(128, 128, 128, 100))  # Increased opacity
                self.rubber_band_center.setWidth(1)  # Increased width
                self.rubber_band_center.setLineStyle(Qt.DashLine)
                
                self.rubber_band1 = QgsRubberBand(self.canvas, QgsWkbTypes.LineGeometry)
                self.rubber_band1.setColor(QColor(255, 0, 0, 255))  # Full opacity
                self.rubber_band1.setWidth(3)  # Increased width
                
                self.rubber_band2 = QgsRubberBand(self.canvas, QgsWkbTypes.LineGeometry)
                self.rubber_band2.setColor(QColor(0, 0, 255, 255))  # Full opacity
                self.rubber_band2.setWidth(3)  # Increased width
                
                # Add initial point to make them visible
                self.rubber_band_center.addPoint(self.start_point)
                self.rubber_band1.addPoint(self.start_point)
                self.rubber_band2.addPoint(self.start_point)
                
                # Ensure they are on top
                self.rubber_band_center.updatePosition()
                self.rubber_band1.updatePosition()
                self.rubber_band2.updatePosition()
            else:
                # Second click: complete drawing
                self.end_point = point
                self.finalize_drawing()
                
    def canvasMoveEvent(self, event):
        """Update visualization while moving mouse"""
        if self.is_drawing and self.start_point:
            current_point = self.toMapCoordinates(event.pos())
            
            # Calculate parallel lines
            line1, line2 = self.calculate_parallel_lines(
                self.start_point, current_point, self.offset_distance
            )
            
            # Update rubber bands
            if self.rubber_band_center:
                self.rubber_band_center.reset(QgsWkbTypes.LineGeometry)
                self.rubber_band_center.addPoint(self.start_point)
                self.rubber_band_center.addPoint(current_point)
                self.rubber_band_center.show()  # Ensure visibility
            
            if self.rubber_band1:
                self.rubber_band1.reset(QgsWkbTypes.LineGeometry)
                for point in line1:
                    self.rubber_band1.addPoint(point)
                self.rubber_band1.show()  # Ensure visibility
                    
            if self.rubber_band2:
                self.rubber_band2.reset(QgsWkbTypes.LineGeometry)
                for point in line2:
                    self.rubber_band2.addPoint(point)
                self.rubber_band2.show()  # Ensure visibility
            
            # Force canvas refresh
            self.canvas.refresh()
                    
    def calculate_parallel_lines(self, start, end, offset):
        """Calculate two parallel lines equidistant from the center line"""
        # Calculate line angle
        dx = end.x() - start.x()
        dy = end.y() - start.y()
        angle = math.atan2(dy, dx)
        
        # Perpendicular angle
        perp_angle = angle + math.pi / 2
        
        # Calculate offset
        offset_x = offset * math.cos(perp_angle)
        offset_y = offset * math.sin(perp_angle)
        
        # Line 1 (A-A' - upper/left)
        line1_start = QgsPointXY(start.x() + offset_x, start.y() + offset_y)
        line1_end = QgsPointXY(end.x() + offset_x, end.y() + offset_y)
        
        # Line 2 (B-B' - lower/right)
        line2_start = QgsPointXY(start.x() - offset_x, start.y() - offset_y)
        line2_end = QgsPointXY(end.x() - offset_x, end.y() - offset_y)
        
        return [line1_start, line1_end], [line2_start, line2_end]
    
    def finalize_drawing(self):
        """Complete drawing and call callback"""
        if self.start_point and self.end_point:
            # Calculate final lines
            line1, line2 = self.calculate_parallel_lines(
                self.start_point, self.end_point, self.offset_distance
            )
            
            # Pass lines to callback
            self.callback(line1, line2, [self.start_point, self.end_point])
            
        # Clean up
        self.cleanup()
        
    def cleanup(self):
        """Clean up resources"""
        for rubber_band in [self.rubber_band_center, self.rubber_band1, self.rubber_band2]:
            if rubber_band:
                self.canvas.scene().removeItem(rubber_band)
        self.rubber_band_center = None
        self.rubber_band1 = None
        self.rubber_band2 = None
        self.start_point = None
        self.end_point = None
        self.is_drawing = False
        
    def deactivate(self):
        """Called when tool is deactivated"""
        self.cleanup()
        super().deactivate()
        
    def keyPressEvent(self, event):
        """Handle ESC to cancel"""
        if event.key() == Qt.Key_Escape:
            self.cleanup()

class ProfileViewerDialog(QtWidgets.QDialog):
    """Main dialog for displaying profiles"""
    
    def __init__(self, iface, parent=None):
        super().__init__(parent)
        self.iface = iface
        self.setWindowTitle("Dual Profile Viewer - Archaeological Analysis")
        self.setMinimumSize(1000, 700)
        
        self.setup_ui()
        
    def setup_ui(self):
        """Create the interface"""
        layout = QtWidgets.QVBoxLayout(self)
        
        # Toolbar
        toolbar = QtWidgets.QHBoxLayout()
        
        # DEM selection
        toolbar.addWidget(QtWidgets.QLabel("Primary DEM:"))
        self.combo_dem = QtWidgets.QComboBox()
        self.load_raster_layers(self.combo_dem)
        toolbar.addWidget(self.combo_dem)
        
        # Multi-DEM comparison
        self.check_multi_dem = QtWidgets.QCheckBox("Compare Multiple DEMs")
        self.check_multi_dem.toggled.connect(self.on_multi_dem_toggled)
        toolbar.addWidget(self.check_multi_dem)
        
        # Add multi-DEM widget (initially hidden)
        self.multi_dem_widget = MultiDEMWidget()
        self.multi_dem_widget.setVisible(False)
        self.multi_dem_widget.selection_changed.connect(self.on_dem_selection_changed)
        layout.addWidget(self.multi_dem_widget)
        
        # Offset distance
        toolbar.addWidget(QtWidgets.QLabel("Offset:"))
        self.spin_distance = QtWidgets.QDoubleSpinBox()
        self.spin_distance.setMinimum(0.05)
        self.spin_distance.setMaximum(500.0)
        self.spin_distance.setValue(1.0)
        self.spin_distance.setSingleStep(0.05)
        self.spin_distance.setDecimals(2)
        self.spin_distance.setSuffix(" m")
        toolbar.addWidget(self.spin_distance)
        
        # Samples
        toolbar.addWidget(QtWidgets.QLabel("Samples:"))
        self.spin_samples = QtWidgets.QSpinBox()
        self.spin_samples.setMinimum(10)
        self.spin_samples.setMaximum(2000)
        self.spin_samples.setValue(200)
        toolbar.addWidget(self.spin_samples)
        
        toolbar.addStretch()
        
        # Action buttons
        self.btn_draw = QtWidgets.QPushButton("🎯 Draw Sections")
        self.btn_draw.clicked.connect(self.start_drawing)
        toolbar.addWidget(self.btn_draw)
        
        self.btn_create_layer = QtWidgets.QPushButton("📍 Create Section Layer")
        self.btn_create_layer.clicked.connect(self.create_section_layer)
        self.btn_create_layer.setEnabled(False)
        toolbar.addWidget(self.btn_create_layer)
        
        self.btn_export_vector = QtWidgets.QPushButton("📐 Export Vector Profile")
        self.btn_export_vector.clicked.connect(self.export_vector_profile)
        self.btn_export_vector.setEnabled(False)
        toolbar.addWidget(self.btn_export_vector)
        
        self.btn_export_data = QtWidgets.QPushButton("💾 Export CSV")
        self.btn_export_data.clicked.connect(self.export_data)
        self.btn_export_data.setEnabled(False)
        toolbar.addWidget(self.btn_export_data)
        
        toolbar.addWidget(QtWidgets.QWidget())  # Spacer
        
        self.btn_help = QtWidgets.QPushButton("❓ Help")
        self.btn_help.clicked.connect(self.show_help)
        self.btn_help.setToolTip("Open user manual")
        toolbar.addWidget(self.btn_help)
        
        layout.addLayout(toolbar)
        
        # Graph area
        info_label = QtWidgets.QLabel(
            "📊 Dual Profile Viewer for Archaeological Analysis\n\n"
            "✅ Graphs will open automatically in your browser\n"
            "✅ Export profiles as georeferenced vectors\n"
            "✅ Multi-DEM comparison support\n\n"
            "Click 'Draw Sections' to begin analysis"
        )
        info_label.setAlignment(qt_compat.AlignCenter)
        info_label.setMinimumHeight(400)
        info_label.setStyleSheet("""
            QLabel { 
                background: qlineargradient(x1:0, y1:0, x2:0, y2:1,
                    stop:0 #f0f4f8, stop:1 #d9e2ec);
                border: 2px solid #627d98;
                border-radius: 8px;
                font-size: 14px;
                padding: 20px;
                color: #102a43;
            }
        """)
        layout.addWidget(info_label)
        
        # Info text area
        self.info_text = QtWidgets.QTextEdit()
        self.info_text.setMaximumHeight(120)
        self.info_text.setReadOnly(True)
        layout.addWidget(self.info_text)
        
        # Variables
        self.profile_data = None
        self.profile_data_list = []  # For 3D viewer
        self.map_tool = None
        self.section_count = 0
        self.current_figure = None
        
    def load_raster_layers(self, combo=None):
        """Load available raster layers"""
        if combo is None:
            combo = self.combo_dem
            
        combo.clear()
        
        for layer in QgsProject.instance().mapLayers().values():
            if isinstance(layer, QgsRasterLayer) and layer.bandCount() > 0:
                combo.addItem(layer.name(), layer.id())
                
        if combo.count() == 0:
            combo.addItem("(No DEM/DTM found)", None)
            
    def on_multi_dem_toggled(self, checked):
        """Toggle visibility of multi-DEM selector"""
        self.multi_dem_widget.setVisible(checked)
        if checked:
            self.multi_dem_widget.refresh_dem_list()
            
    def on_dem_selection_changed(self, selected_layers):
        """Handle changes in DEM selection"""
        # Store selected DEMs for comparison
        self.selected_dem_layers = selected_layers
            
    def start_drawing(self):
        """Activate drawing tool"""
        if self.combo_dem.currentData() is None:
            QtWidgets.QMessageBox.warning(
                self, 
                "Warning", 
                "Please select a DEM/DTM layer first!"
            )
            return
            
        # Set distance in tool
        canvas = self.iface.mapCanvas()
        
        # Deactivate previous tool if exists
        if self.map_tool:
            canvas.unsetMapTool(self.map_tool)
            
        # Create new tool
        self.map_tool = DualProfileTool(canvas, self.on_lines_drawn)
        self.map_tool.offset_distance = self.spin_distance.value()
        
        # Activate tool
        canvas.setMapTool(self.map_tool)
        
        self.info_text.setText("Click two points on the map to draw parallel sections...\n"
                              "The sections will be created parallel to and equidistant from your drawn line.")
        
    def on_lines_drawn(self, line1, line2, center_line):
        """Called when lines have been drawn"""
        # Deactivate tool
        self.iface.mapCanvas().unsetMapTool(self.map_tool)
        
        # Extract profiles
        self.extract_profiles(line1, line2, center_line)
        
    def extract_profiles(self, line1, line2, center_line):
        """Extract elevation profiles from lines"""
        dem_layer_id = self.combo_dem.currentData()
        if not dem_layer_id:
            return
            
        dem_layer = QgsProject.instance().mapLayer(dem_layer_id)
        if not dem_layer:
            return
            
        num_samples = self.spin_samples.value()
        
        # Extract profiles
        profile1 = self.sample_raster_along_line(
            dem_layer, line1[0], line1[1], num_samples
        )
        
        profile2 = self.sample_raster_along_line(
            dem_layer, line2[0], line2[1], num_samples
        )
        
        # Debug primary DEM
        valid_count1 = np.sum(~np.isnan(profile1['elevations']))
        print(f"DEM1 ({dem_layer.name()}): {valid_count1}/{num_samples} valid points")
        
        # Check for multiple DEMs
        additional_profiles = {}
        
        if self.check_multi_dem.isChecked() and hasattr(self, 'selected_dem_layers'):
            for i, layer_id in enumerate(self.selected_dem_layers):
                if layer_id == dem_layer.id():
                    continue  # Skip primary DEM
                    
                dem_layer_comp = QgsProject.instance().mapLayer(layer_id)
                if dem_layer_comp:
                    dem_name = dem_layer_comp.name()
                    
                    # Check if DEMs have same CRS
                    if dem_layer.crs() != dem_layer_comp.crs():
                        QtWidgets.QMessageBox.warning(
                            self,
                            "CRS Mismatch",
                            f"Warning: DEMs have different coordinate systems!\n"
                            f"DEM1: {dem_layer.crs().authid()}\n"
                            f"DEM2: {dem2_layer.crs().authid()}\n\n"
                            f"Results may be incorrect. Consider reprojecting one DEM."
                        )
                    
                    profile1_dem2 = self.sample_raster_along_line(
                        dem2_layer, line1[0], line1[1], num_samples
                    )
                    profile2_dem2 = self.sample_raster_along_line(
                        dem2_layer, line2[0], line2[1], num_samples
                    )
                    
                    # Debug secondary DEM
                    valid_count2_a = np.sum(~np.isnan(profile1_dem2['elevations']))
                    valid_count2_b = np.sum(~np.isnan(profile2_dem2['elevations']))
                    print(f"DEM2 ({dem2_name}): A-A' {valid_count2_a}/{num_samples}, B-B' {valid_count2_b}/{num_samples} valid points")
                    
                    if valid_count2_a == 0 and valid_count2_b == 0:
                        QtWidgets.QMessageBox.warning(
                            self,
                            "No Data in Second DEM",
                            f"The second DEM '{dem2_name}' has no valid elevation data\n"
                            f"in the selected area.\n\n"
                            f"Check that:\n"
                            f"1. The DEM covers this area\n"
                            f"2. The DEM has valid elevation values\n"
                            f"3. Both DEMs use the same coordinate system"
                        )
                        # Reset to None to avoid plotting empty data
                        profile1_dem2 = None
                        profile2_dem2 = None
                        dem2_name = None
        
        # Save data
        self.profile_data = {
            'profile1': profile1,
            'profile2': profile2,
            'profile1_dem2': profile1_dem2,
            'profile2_dem2': profile2_dem2,
            'line1': line1,
            'line2': line2,
            'center_line': center_line,
            'offset': self.spin_distance.value(),
            'dem1_name': dem_layer.name(),
            'dem2_name': dem2_name
        }
        
        # Visualize
        self.plot_profiles(profile1, profile2, profile1_dem2, profile2_dem2)
        
        # Enable export buttons
        self.btn_export_data.setEnabled(True)
        self.btn_export_vector.setEnabled(True)
        self.btn_create_layer.setEnabled(True)
        
    def sample_raster_along_line(self, raster_layer, start_point, end_point, num_samples):
        """Sample raster along a line with better NoData handling"""
        # Create interpolated points along line
        x_coords = np.linspace(start_point.x(), end_point.x(), num_samples)
        y_coords = np.linspace(start_point.y(), end_point.y(), num_samples)
        
        # Calculate progressive distances
        distances = np.zeros(num_samples)
        elevations = np.zeros(num_samples)
        
        # Get NoData value for this raster
        provider = raster_layer.dataProvider()
        nodata_value = provider.sourceNoDataValue(1)  # Band 1
        
        # Check if points are within raster extent
        extent = raster_layer.extent()
        points_in_extent = 0
        
        for i in range(num_samples):
            if i > 0:
                dx = x_coords[i] - x_coords[i-1]
                dy = y_coords[i] - y_coords[i-1]
                distances[i] = distances[i-1] + math.sqrt(dx*dx + dy*dy)
            
            # Sample elevation
            sample_point = QgsPointXY(x_coords[i], y_coords[i])
            
            # Check if point is within raster extent
            if not extent.contains(sample_point):
                elevations[i] = np.nan
                continue
            
            points_in_extent += 1
            
            # Sample the raster value
            value, success = provider.sample(sample_point, 1)
            
            if success and value != nodata_value and not math.isnan(value):
                elevations[i] = value
            else:
                elevations[i] = np.nan
        
        # Debug output
        valid_count = np.sum(~np.isnan(elevations))
        if valid_count == 0:
            print(f"WARNING: No valid elevation data extracted from {raster_layer.name()}")
            print(f"  Points in extent: {points_in_extent}/{num_samples}")
            print(f"  NoData value: {nodata_value}")
            if points_in_extent > 0:
                print(f"  All values were NoData or failed to sample")
                
        return {
            'distances': distances,
            'elevations': elevations,
            'x_coords': x_coords,
            'y_coords': y_coords
        }
        
    def plot_profiles(self, profile1, profile2, profile1_dem2=None, profile2_dem2=None):
        """Display profiles with multi-DEM support"""
        has_multi_dem = profile1_dem2 is not None and profile2_dem2 is not None
        
        if PLOTLY_AVAILABLE:
            # Create Plotly figure
            fig = make_subplots(
                rows=2, cols=1,
                subplot_titles=(
                    f'Parallel Profiles - Offset: {self.spin_distance.value():.2f}m',
                    "Elevation Difference (A-A' minus B-B')"
                ),
                vertical_spacing=0.15,
                row_heights=[0.6, 0.4]
            )
            
            # Add primary DEM profiles
            fig.add_trace(
                go.Scatter(
                    x=profile1['distances'], 
                    y=profile1['elevations'],
                    mode='lines',
                    name=f"A-A' ({self.profile_data['dem1_name']})",
                    line=dict(color='red', width=2),
                    hovertemplate='Distance: %{x:.2f}m<br>Elevation: %{y:.3f}m<extra></extra>'
                ),
                row=1, col=1
            )
            
            fig.add_trace(
                go.Scatter(
                    x=profile2['distances'], 
                    y=profile2['elevations'],
                    mode='lines',
                    name=f"B-B' ({self.profile_data['dem1_name']})",
                    line=dict(color='blue', width=2),
                    hovertemplate='Distance: %{x:.2f}m<br>Elevation: %{y:.3f}m<extra></extra>'
                ),
                row=1, col=1
            )
            
            # Add fill between primary profiles
            fig.add_trace(
                go.Scatter(
                    x=np.concatenate([profile1['distances'], profile2['distances'][::-1]]),
                    y=np.concatenate([profile1['elevations'], profile2['elevations'][::-1]]),
                    fill='toself',
                    fillcolor='rgba(128, 128, 128, 0.2)',
                    line=dict(color='rgba(255,255,255,0)'),
                    name='Volume/Difference',
                    hoverinfo='skip',
                    showlegend=True
                ),
                row=1, col=1
            )
            
            # Add secondary DEM profiles if available
            if has_multi_dem:
                fig.add_trace(
                    go.Scatter(
                        x=profile1_dem2['distances'], 
                        y=profile1_dem2['elevations'],
                        mode='lines',
                        name=f"A-A' ({self.profile_data['dem2_name']})",
                        line=dict(color='red', width=1.5, dash='dash'),
                        hovertemplate='Distance: %{x:.2f}m<br>Elevation: %{y:.3f}m<extra></extra>'
                    ),
                    row=1, col=1
                )
                
                fig.add_trace(
                    go.Scatter(
                        x=profile2_dem2['distances'], 
                        y=profile2_dem2['elevations'],
                        mode='lines',
                        name=f"B-B' ({self.profile_data['dem2_name']})",
                        line=dict(color='blue', width=1.5, dash='dash'),
                        hovertemplate='Distance: %{x:.2f}m<br>Elevation: %{y:.3f}m<extra></extra>'
                    ),
                    row=1, col=1
                )
            
            # Calculate and plot difference for primary DEM
            diff1 = profile1['elevations'] - profile2['elevations']
            
            fig.add_trace(
                go.Scatter(
                    x=profile1['distances'], 
                    y=diff1,
                    mode='lines',
                    name=f'Difference ({self.profile_data["dem1_name"]})',
                    line=dict(color='green', width=2),
                    hovertemplate='Distance: %{x:.2f}m<br>Difference: %{y:.3f}m<extra></extra>'
                ),
                row=2, col=1
            )
            
            # Add fills for positive/negative differences
            fig.add_trace(
                go.Scatter(
                    x=profile1['distances'],
                    y=np.where(diff1 >= 0, diff1, 0),
                    fill='tozeroy',
                    fillcolor='rgba(0, 255, 0, 0.3)',
                    line=dict(color='rgba(255,255,255,0)'),
                    name='A higher',
                    hoverinfo='skip',
                    showlegend=True
                ),
                row=2, col=1
            )
            
            fig.add_trace(
                go.Scatter(
                    x=profile1['distances'],
                    y=np.where(diff1 < 0, diff1, 0),
                    fill='tozeroy',
                    fillcolor='rgba(255, 0, 0, 0.3)',
                    line=dict(color='rgba(255,255,255,0)'),
                    name='B higher',
                    hoverinfo='skip',
                    showlegend=True
                ),
                row=2, col=1
            )
            
            # Add secondary DEM difference if available
            if has_multi_dem:
                diff2 = profile1_dem2['elevations'] - profile2_dem2['elevations']
                fig.add_trace(
                    go.Scatter(
                        x=profile1_dem2['distances'], 
                        y=diff2,
                        mode='lines',
                        name=f'Difference ({self.profile_data["dem2_name"]})',
                        line=dict(color='orange', width=2, dash='dash'),
                        hovertemplate='Distance: %{x:.2f}m<br>Difference: %{y:.3f}m<extra></extra>'
                    ),
                    row=2, col=1
                )
            
            # Add zero line
            fig.add_trace(
                go.Scatter(
                    x=[profile1['distances'][0], profile1['distances'][-1]],
                    y=[0, 0],
                    mode='lines',
                    line=dict(color='black', width=1, dash='dash'),
                    name='Zero line',
                    showlegend=False,
                    hoverinfo='skip'
                ),
                row=2, col=1
            )
            
            # Update layout
            fig.update_xaxes(title_text="Distance (m)", row=1, col=1, showgrid=True)
            fig.update_xaxes(title_text="Distance (m)", row=2, col=1, showgrid=True)
            fig.update_yaxes(title_text="Elevation (m)", row=1, col=1, showgrid=True)
            fig.update_yaxes(title_text="Elevation Difference (m)", row=2, col=1, showgrid=True)
            
            fig.update_layout(
                height=600,
                showlegend=True,
                hovermode='x unified',
                template='plotly_white'
            )
            
            # Store figure
            self.current_figure = fig
            
            # Open in browser
            self.open_graph_in_browser(fig)
            
        # Calculate statistics
        if has_multi_dem:
            self.calculate_statistics_multi(profile1, profile2, profile1_dem2, profile2_dem2)
        else:
            self.calculate_statistics(profile1, profile2)
    
    def open_graph_in_browser(self, fig):
        """Open graph in browser"""
        temp_html = os.path.join(tempfile.gettempdir(), f"profile_{datetime.now().strftime('%Y%m%d_%H%M%S')}.html")
        fig.write_html(temp_html, auto_open=True)
        
        self.info_text.setText(
            f"📊 INTERACTIVE GRAPH OPENED IN BROWSER\n"
            f"File: {temp_html}\n\n"
            f"Use the browser's tools to interact with the graph."
        )
    
    def calculate_statistics(self, profile1, profile2):
        """Calculate and show statistics"""
        diff = profile1['elevations'] - profile2['elevations']
        valid_diff = diff[~np.isnan(diff)]
        
        if len(valid_diff) > 0:
            stats_text = f"""📊 PROFILE STATISTICS:
━━━━━━━━━━━━━━━━━━━━━
Profile A-A' (Upper):
  • Min: {np.nanmin(profile1['elevations']):.3f} m
  • Max: {np.nanmax(profile1['elevations']):.3f} m
  • Mean: {np.nanmean(profile1['elevations']):.3f} m

Profile B-B' (Lower):
  • Min: {np.nanmin(profile2['elevations']):.3f} m
  • Max: {np.nanmax(profile2['elevations']):.3f} m
  • Mean: {np.nanmean(profile2['elevations']):.3f} m

Difference (A-B):
  • Min: {np.min(valid_diff):.3f} m
  • Max: {np.max(valid_diff):.3f} m
  • Mean: {np.mean(valid_diff):.3f} m

Section length: {profile1['distances'][-1]:.2f} m
Offset: {self.spin_distance.value():.2f} m
            """
        else:
            stats_text = "No valid data extracted"
            
        self.info_text.setText(stats_text)
    
    def calculate_statistics_multi(self, profile1, profile2, profile1_dem2, profile2_dem2):
        """Calculate and show statistics for multiple DEMs"""
        diff1 = profile1['elevations'] - profile2['elevations']
        diff2 = profile1_dem2['elevations'] - profile2_dem2['elevations']
        
        # Remove NaN for statistics
        valid_diff1 = diff1[~np.isnan(diff1)]
        valid_diff2 = diff2[~np.isnan(diff2)]
        
        # DEM differences
        dem_diff_a = profile1['elevations'] - profile1_dem2['elevations']
        dem_diff_b = profile2['elevations'] - profile2_dem2['elevations']
        valid_dem_diff_a = dem_diff_a[~np.isnan(dem_diff_a)]
        valid_dem_diff_b = dem_diff_b[~np.isnan(dem_diff_b)]
        
        stats_text = f"""📊 MULTI-DEM COMPARISON:
━━━━━━━━━━━━━━━━━━━━━
{self.profile_data['dem1_name']}:
  A-A': Min={np.nanmin(profile1['elevations']):.3f}, Max={np.nanmax(profile1['elevations']):.3f}, Mean={np.nanmean(profile1['elevations']):.3f}
  B-B': Min={np.nanmin(profile2['elevations']):.3f}, Max={np.nanmax(profile2['elevations']):.3f}, Mean={np.nanmean(profile2['elevations']):.3f}
  Diff: Min={np.min(valid_diff1) if len(valid_diff1) > 0 else 0:.3f}, Max={np.max(valid_diff1) if len(valid_diff1) > 0 else 0:.3f}, Mean={np.mean(valid_diff1) if len(valid_diff1) > 0 else 0:.3f}

{self.profile_data['dem2_name']}:
  A-A': Min={np.nanmin(profile1_dem2['elevations']):.3f}, Max={np.nanmax(profile1_dem2['elevations']):.3f}, Mean={np.nanmean(profile1_dem2['elevations']):.3f}
  B-B': Min={np.nanmin(profile2_dem2['elevations']):.3f}, Max={np.nanmax(profile2_dem2['elevations']):.3f}, Mean={np.nanmean(profile2_dem2['elevations']):.3f}
  Diff: Min={np.min(valid_diff2) if len(valid_diff2) > 0 else 0:.3f}, Max={np.max(valid_diff2) if len(valid_diff2) > 0 else 0:.3f}, Mean={np.mean(valid_diff2) if len(valid_diff2) > 0 else 0:.3f}

DEM Comparison:
  A-A' difference: Mean={np.mean(valid_dem_diff_a) if len(valid_dem_diff_a) > 0 else 0:.3f}m, StdDev={np.std(valid_dem_diff_a) if len(valid_dem_diff_a) > 0 else 0:.3f}m
  B-B' difference: Mean={np.mean(valid_dem_diff_b) if len(valid_dem_diff_b) > 0 else 0:.3f}m, StdDev={np.std(valid_dem_diff_b) if len(valid_dem_diff_b) > 0 else 0:.3f}m

Section: {profile1['distances'][-1]:.2f}m, Points: {len(profile1['distances'])}, Offset: {self.spin_distance.value():.2f}m
        """
        
        self.info_text.setText(stats_text)
        
    def create_section_layer(self):
        """Create a vector layer with section lines"""
        if not self.profile_data:
            return
            
        self.section_count += 1
        
        # Create or get existing layer
        layer_name = "Profile_Sections"
        layers = QgsProject.instance().mapLayersByName(layer_name)
        
        if layers:
            section_layer = layers[0]
        else:
            # Create new layer
            crs = QgsProject.instance().crs()
            section_layer = QgsVectorLayer(
                f"LineString?crs={crs.authid()}&field=id:integer&field=label:string&field=type:string",
                layer_name,
                "memory"
            )
            
            QgsProject.instance().addMapLayer(section_layer)
        
        # Add features
        section_layer.startEditing()
        
        # Add line A-A'
        feature1 = QgsFeature()
        feature1.setGeometry(QgsGeometry.fromPolylineXY(self.profile_data['line1']))
        feature1.setAttributes([
            self.section_count * 2 - 1,
            f"A{self.section_count}-A'{self.section_count}",
            "Upper"
        ])
        section_layer.addFeature(feature1)
        
        # Add line B-B'
        feature2 = QgsFeature()
        feature2.setGeometry(QgsGeometry.fromPolylineXY(self.profile_data['line2']))
        feature2.setAttributes([
            self.section_count * 2,
            f"B{self.section_count}-B'{self.section_count}",
            "Lower"
        ])
        section_layer.addFeature(feature2)
        
        section_layer.commitChanges()
        
        # Apply symbology
        self.setup_layer_symbology(section_layer)
        
        self.iface.mapCanvas().refresh()
        
        QtWidgets.QMessageBox.information(
            self,
            "Success",
            f"✅ Section lines created"
        )
        
    def setup_layer_symbology(self, layer):
        """Set up symbology"""
        categories = []
        
        # Red for Upper
        symbol_upper = QgsSymbol.defaultSymbol(layer.geometryType())
        symbol_upper.setColor(QColor(255, 0, 0))
        symbol_upper.setWidth(1.5)
        cat_upper = QgsRendererCategory('Upper', symbol_upper, "A-A'")
        categories.append(cat_upper)
        
        # Blue for Lower
        symbol_lower = QgsSymbol.defaultSymbol(layer.geometryType())
        symbol_lower.setColor(QColor(0, 0, 255))
        symbol_lower.setWidth(1.5)
        cat_lower = QgsRendererCategory('Lower', symbol_lower, "B-B'")
        categories.append(cat_lower)
        
        renderer = QgsCategorizedSymbolRenderer('type', categories)
        layer.setRenderer(renderer)
        layer.triggerRepaint()
    
    def export_vector_profile(self):
        """Export profile as georeferenced vector"""
        if not self.profile_data:
            QtWidgets.QMessageBox.warning(self, "Warning", 
                                         "No profile data to export!")
            return
        
        # Show options dialog
        dialog = VectorExportDialog(self)
        if qt_compat.dialog_exec(dialog):
            options = dialog.get_options()
            
            # Check if we should include secondary DEM profiles
            include_dem2 = False
            if self.profile_data.get('profile1_dem2') is not None:
                # Check if secondary DEM has valid data
                valid_count = np.sum(~np.isnan(self.profile_data['profile1_dem2']['elevations']))
                if valid_count > 0:
                    include_dem2 = True
                    
                    # Ask user if they want to export both DEMs
                    reply = QtWidgets.QMessageBox.question(
                        self,
                        "Export Multiple DEMs",
                        "You have compared two DEMs.\n\n"
                        "Do you want to export profiles from both DEMs?\n"
                        "Yes = Export both\n"
                        "No = Export primary DEM only",
                        QtWidgets.QMessageBox.Yes | QtWidgets.QMessageBox.No
                    )
                    
                    if reply == QtWidgets.QMessageBox.No:
                        include_dem2 = False
            
            # Get output file
            filename, _ = QtWidgets.QFileDialog.getSaveFileName(
                self,
                "Save profile as vector",
                f"profile_vector_{datetime.now().strftime('%Y%m%d_%H%M%S')}.gpkg",
                "GeoPackage (*.gpkg);;Shapefile (*.shp)"
            )
            
            if filename:
                try:
                    # Handle 3D export separately
                    if options['type'] == '3d':
                        layer = ProfileExporter.export_profiles_as_3d_vector(
                            self.profile_data,
                            filename,
                            add_to_map=options['add_to_map']
                        )
                    else:
                        # For 2D export, handle multiple DEMs
                        if include_dem2:
                            # Export primary DEM
                            base_name = os.path.splitext(filename)[0]
                            ext = os.path.splitext(filename)[1]
                            
                            # Create filenames for each DEM
                            filename1 = f"{base_name}_DEM1{ext}"
                            filename2 = f"{base_name}_DEM2{ext}"
                            
                            # Export primary DEM profiles
                            profile_data_dem1 = {
                                'profile1': self.profile_data['profile1'],
                                'profile2': self.profile_data['profile2'],
                                'line1': self.profile_data['line1'],
                                'line2': self.profile_data['line2']
                            }
                            
                            layer1 = ProfileExporter.export_profile_as_vector(
                                profile_data_dem1,
                                filename1,
                                export_type=options['type'],
                                scale_factor=options['scale'],
                                vertical_exaggeration=options['vertical_exag'],
                                baseline_offset=options['baseline_offset'],
                                add_to_map=options['add_to_map']
                            )
                            
                            # Export secondary DEM profiles
                            profile_data_dem2 = {
                                'profile1': self.profile_data['profile1_dem2'],
                                'profile2': self.profile_data['profile2_dem2'],
                                'line1': self.profile_data['line1'],
                                'line2': self.profile_data['line2']
                            }
                            
                            layer2 = ProfileExporter.export_profile_as_vector(
                                profile_data_dem2,
                                filename2,
                                export_type=options['type'],
                                scale_factor=options['scale'],
                                vertical_exaggeration=options['vertical_exag'],
                                baseline_offset=options['baseline_offset'],
                                add_to_map=options['add_to_map']
                            )
                            
                            QtWidgets.QMessageBox.information(
                                self, "Success",
                                f"✅ Profiles exported from both DEMs:\n"
                                f"DEM1: {filename1}\n"
                                f"DEM2: {filename2}\n\n"
                                f"The profile shapes have been preserved as georeferenced features."
                            )
                        else:
                            # Export single DEM
                            layer = ProfileExporter.export_profile_as_vector(
                                self.profile_data,
                                filename,
                                export_type=options['type'],
                                scale_factor=options['scale'],
                                vertical_exaggeration=options['vertical_exag'],
                                baseline_offset=options['baseline_offset'],
                                add_to_map=options['add_to_map']
                            )
                            
                            QtWidgets.QMessageBox.information(
                                self, "Success",
                                f"✅ Profile exported as vector:\n{filename}\n\n"
                                f"The profile shape has been preserved as a georeferenced feature."
                            )
                    
                except Exception as e:
                    QtWidgets.QMessageBox.critical(
                        self, "Error", f"Export failed:\n{str(e)}"
                    )
    
    def export_data(self):
        """Export profile data to CSV"""
        if not self.profile_data:
            return
            
        filename, _ = QtWidgets.QFileDialog.getSaveFileName(
            self,
            "Save profiles",
            f"profiles_{datetime.now().strftime('%Y%m%d_%H%M%S')}.csv",
            "CSV Files (*.csv)"
        )
        
        if filename:
            try:
                import csv
                
                with open(filename, 'w', newline='') as csvfile:
                    writer = csv.writer(csvfile)
                    
                    # Header
                    header = [
                        'Distance_m', 'X', 'Y',
                        'Elevation_A_m', 'Elevation_B_m',
                        'Difference_m', 'Offset_m'
                    ]
                    writer.writerow(header)
                    
                    # Data
                    p1 = self.profile_data['profile1']
                    p2 = self.profile_data['profile2']
                    
                    for i in range(len(p1['distances'])):
                        row = [
                            p1['distances'][i],
                            p1['x_coords'][i],
                            p1['y_coords'][i],
                            p1['elevations'][i],
                            p2['elevations'][i],
                            p1['elevations'][i] - p2['elevations'][i],
                            self.profile_data['offset']
                        ]
                        writer.writerow(row)
                        
                QtWidgets.QMessageBox.information(
                    self, "Success",
                    f"Data exported to:\n{filename}"
                )
                
            except Exception as e:
                QtWidgets.QMessageBox.critical(
                    self, "Error",
                    f"Export error:\n{str(e)}"
                )
    
    def show_help(self):
        """Open the user manual"""
        import os
        import webbrowser
        from pathlib import Path
        
        # Get the plugin directory
        plugin_dir = os.path.dirname(__file__)
        manual_path = os.path.join(plugin_dir, "USER_MANUAL.md")
        
        # Check if manual exists
        if os.path.exists(manual_path):
            # Try to open with system default markdown viewer
            # On most systems, this will open in the default browser or markdown editor
            try:
                # Convert to file URL
                file_url = Path(manual_path).as_uri()
                webbrowser.open(file_url)
            except Exception:
                # If that fails, try to open directly
                try:
                    import subprocess
                    if os.name == 'nt':  # Windows
                        os.startfile(manual_path)
                    elif os.name == 'posix':  # macOS and Linux
                        subprocess.call(['open', manual_path])
                except Exception as e:
                    # Last resort: show path to user
                    QtWidgets.QMessageBox.information(
                        self,
                        "User Manual",
                        f"User manual is located at:\n\n{manual_path}\n\n"
                        f"Please open it with your preferred markdown viewer."
                    )
        else:
            # Manual not found, offer to open online version
            reply = QtWidgets.QMessageBox.question(
                self,
                "User Manual",
                "Local manual not found.\n\n"
                "Would you like to open the online documentation?",
                QtWidgets.QMessageBox.Yes | QtWidgets.QMessageBox.No
            )
            
            if reply == QtWidgets.QMessageBox.Yes:
                webbrowser.open("https://github.com/enzococca/dual-profile-viewer/wiki")
    
    def open_3d_viewer(self):
        """Open the 3D viewer with current profile data"""
        if not self.profile_data:
            QtWidgets.QMessageBox.warning(self, "No Data", 
                                         "Please create profiles first before opening 3D viewer")
            return
        
        # Import the 3D viewer
        from .advanced_3d_viewer import Advanced3DViewer
        
        # Create viewer if not exists
        if not hasattr(self, 'viewer_3d') or self.viewer_3d is None:
            self.viewer_3d = Advanced3DViewer(self)
        
        # Prepare profile data list with proper format
        self.prepare_profile_data_for_3d()
        
        # Load profiles into viewer
        if self.profile_data_list:
            self.viewer_3d.load_profiles(self.profile_data_list)
        
        # Show the viewer
        self.viewer_3d.show()
    
    def prepare_profile_data_for_3d(self):
        """Prepare profile data in format needed for 3D viewer"""
        self.profile_data_list = []
        
        if not self.profile_data:
            return
        
        # Prepare profile 1
        if 'profile1' in self.profile_data:
            p1 = self.profile_data['profile1']
            profile1_data = {
                'name': 'Profile A (Red)',
                'distances': p1['distances'],
                'elevations': p1['elevations'],
                'coordinates': list(zip(p1['x_coords'], p1['y_coords'])) if 'x_coords' in p1 else None
            }
            self.profile_data_list.append(profile1_data)
        
        # Prepare profile 2
        if 'profile2' in self.profile_data:
            p2 = self.profile_data['profile2']
            profile2_data = {
                'name': 'Profile B (Blue)',
                'distances': p2['distances'],
                'elevations': p2['elevations'],
                'coordinates': list(zip(p2['x_coords'], p2['y_coords'])) if 'x_coords' in p2 else None
            }
            self.profile_data_list.append(profile2_data)
    
    def get_profile_data(self):
        """Get profile data for 3D visualization"""
        if not self.profile_data:
            return None
            
        # Return combined profile data suitable for stratigraph
        profile_points = []
        
        # Use profile 1 as the main profile
        p1 = self.profile_data['profile1']
        for i in range(len(p1['distances'])):
            if not np.isnan(p1['elevations'][i]):
                profile_points.append((
                    p1['distances'][i],
                    p1['elevations'][i],
                    0  # Z coordinate for 2D profile
                ))
                
        return profile_points
        
    def get_all_profile_data(self):
        """Get all profile data for cross-section visualization"""
        if not self.profile_data:
            return []
            
        all_profiles = []
        
        # Get both profiles
        for profile_key in ['profile1', 'profile2']:
            if profile_key in self.profile_data:
                profile = self.profile_data[profile_key]
                profile_points = []
                
                for i in range(len(profile['distances'])):
                    if not np.isnan(profile['elevations'][i]):
                        profile_points.append((
                            profile['distances'][i],
                            profile['elevations'][i],
                            0
                        ))
                        
                if profile_points:
                    all_profiles.append(profile_points)
                    
        return all_profiles