# -*- coding: utf-8 -*-

"""
/***************************************************************************
 Geohazard
                               A QGIS plugin
 Plugin with various tools for landslide analysis and management
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                               -------------------
        begin                : 2024-11-29
        copyright            : (C) 2024 by Campus S., Castelli M., Fasciano C., Filipello A.
        email                : andrea.filipello@arpa.piemonte.it
 ***************************************************************************/

/***************************************************************************
 * *
 * This program is free software; you can redistribute it and/or modify  *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation; either version 2 of the License, or     *
 * (at your option) any later version.                                   *
 * *
 * *************************************************************************/
"""

__author__ = 'Campus S., Castelli M., Fasciano C., Filipello A.'
__date__ = '2024-11-29'
__copyright__ = '(C) 2024 by Campus S., Castelli M., Fasciano C., Filipello A.'

# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'

import os
from PyQt5.QtGui import QIcon 
from qgis.core import (
    QgsProcessing,
    QgsProcessingAlgorithm,
    QgsProcessingMultiStepFeedback,
    QgsProcessingParameterRasterLayer,
    QgsProcessingParameterNumber,
    QgsProcessingParameterRasterDestination,
    QgsProcessingLayerPostProcessorInterface,
    QgsColorRampShader,
    QgsProcessingContext
)
import processing

# Variabile globale che mantiene in memoria l'istanza del PostProcessor impedendo
# al Garbage Collector di Python di distruggerla prima del caricamento del layer.
style_processor_instance = None


class VisibilityStyleProcessor(QgsProcessingLayerPostProcessorInterface):
    """
    Post-Processor personalizzato incaricato di caricare il file .qml,
    forzare l'interfaccia grafica su 'Banda singola falso colore' e impostare
    la rampa su modalità discreta a protezione del raster riclassificato.
    """
    def __init__(self, qml_path):
        super().__init__()
        self.qml_path = qml_path
        

    def postProcessLayer(self, layer, context, feedback):
        if layer and os.path.exists(self.qml_path):
            # 1. Carica le proprietà simbologiche registrate nel file .qml
            success, message = layer.loadNamedStyle(self.qml_path)
            
            if success:
                # 2. Forza QGIS a impostare l'albero delle legende e l'interfaccia nativa
                # sulla modalità 'Banda singola falso colore' anziché sulla scala di grigi di default.
                layer.setHtmlProperty("rendererType", "singlebandpseudocolor")
                
                # 3. Interviene sul raster shader impostando le classi come 'Discrete'.
                # Corregge l'errore del .qml originario (INTERPOLATED) non compatibile con classi intere 1-5.
                renderer = layer.renderer()
                if hasattr(renderer, 'shader'):
                    shader_function = renderer.shader().rasterShaderFunction()
                    if isinstance(shader_function, QgsColorRampShader):
                        shader_function.setColorRampType(QgsColorRampShader.Discrete)
                
                # Forza il ridisegno della Canvas di QGIS con i nuovi parametri grafici
                layer.triggerRepaint()
                if feedback:
                    feedback.pushInfo("Stile .qml applicato e interfaccia Falso Colore sbloccata con successo!")
            else:
                if feedback:
                    feedback.pushInfo(f"Errore nel caricamento del file .qml: {message}")
        else:
            if feedback:
                feedback.pushInfo(f"File di stile non trovato al percorso: {self.qml_path}")


class GroundmotionCIndex(QgsProcessingAlgorithm):

    def createInstance(self):
        # Ritorna una nuova istanza della classe, risolvendo l'errore NotImplementedError
        return GroundmotionCIndex()

    def icon(self):
        """
        Should return a QIcon which is used for your provider inside
        the Processing toolbox.
        """
        cmd_folder = os.path.dirname(__file__)
        icon_path = os.path.join(cmd_folder, 'icons', 'Cindex_algorithm.png')
        return QIcon(icon_path)

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterRasterLayer('dtm', 'DTM', defaultValue=None))
        self.addParameter(QgsProcessingParameterNumber('cell_size', 'Cell size', type=QgsProcessingParameterNumber.Integer, defaultValue=5))
        self.addParameter(QgsProcessingParameterNumber('__incidence_angle_090', 'θ - Incidence angle (0-90)', type=QgsProcessingParameterNumber.Double, minValue=0, maxValue=90, defaultValue=39))
        self.addParameter(QgsProcessingParameterNumber('__track_angle__azimut_0360', 'φ - Track angle - Azimut (0-360)', type=QgsProcessingParameterNumber.Double, minValue=0, maxValue=360, defaultValue=None))
        self.addParameter(QgsProcessingParameterRasterDestination('C_index', 'C_Index', createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterRasterDestination('Percentage_c_index', 'Percentage_c_index', createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterRasterDestination('Class_of_visibility', 'Class_of_visibility', createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        feedback = QgsProcessingMultiStepFeedback(11, model_feedback)
        results = {}
        outputs = {}

        # --- ESTRAZIONE PARAMETRI E COSTRUZIONE NOMI DINAMICI ---
        inc_angle = parameters.get('__incidence_angle_090', 39)
        track_angle = parameters.get('__track_angle__azimut_0360', 0)
        
        # Sostituiamo i punti con underscore per pulizia (es: 39.5 -> 39_5)
        inc_str = str(inc_angle).replace('.', '_')
        track_str = str(track_angle).replace('.', '_')
        
        # Suffisso comune per i nomi dei layer in legenda
        suffix = f"_track{track_str}_inc{inc_str}"
        
        # Nomi visualizzati finali in QGIS
        name_c_index = f"C_Index{suffix}"
        name_percentage = f"Percentage_c_index{suffix}"
        name_visibility = f"Class_of_visibility{suffix}"

        # Gestione dei percorsi fisici se l'utente ha scelto "Salva su file"
        c_index_dest = parameters['C_index']
        if c_index_dest and isinstance(c_index_dest, str) and not c_index_dest.startswith('memory:'):
            base, ext = os.path.splitext(c_index_dest)
            if not base.endswith(suffix):
                c_index_dest = f"{base}{suffix}{ext}"

        perc_c_index_dest = parameters['Percentage_c_index']
        if perc_c_index_dest and isinstance(perc_c_index_dest, str) and not perc_c_index_dest.startswith('memory:'):
            base, ext = os.path.splitext(perc_c_index_dest)
            if not base.endswith(suffix):
                perc_c_index_dest = f"{base}{suffix}{ext}"

        class_vis_dest = parameters['Class_of_visibility']
        if class_vis_dest and isinstance(class_vis_dest, str) and not class_vis_dest.startswith('memory:'):
            base, ext = os.path.splitext(class_vis_dest)
            if not base.endswith(suffix):
                class_vis_dest = f"{base}{suffix}{ext}"
        # --------------------------------------------------------

        # Raster_φ - Track angle
        alg_params = {
            'EXTENT': parameters['dtm'],
            'NUMBER': parameters['__track_angle__azimut_0360'],
            'OUTPUT_TYPE': 5,  # Float32
            'PIXEL_SIZE': parameters['cell_size'],
            'TARGET_CRS': parameters['dtm'],
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['Raster_TrackAngle'] = processing.run('native:createconstantrasterlayer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(1)
        if feedback.isCanceled():
            return {}

        # Azimut angle
        alg_params = {
            'BAND_A': 1,
            'BAND_B': None,
            'BAND_C': None,
            'BAND_D': None,
            'BAND_E': None,
            'BAND_F': None,
            'EXTENT_OPT': 0,
            'EXTRA': '',
            'FORMULA': 'A+90',
            'INPUT_A': outputs['Raster_TrackAngle']['OUTPUT'],
            'INPUT_B': None,
            'INPUT_C': None,
            'INPUT_D': None,
            'INPUT_E': None,
            'INPUT_F': None,
            'NO_DATA': None,
            'OPTIONS': '',
            'PROJWIN': None,
            'RTYPE': 5,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['AzimutAngle'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(2)
        if feedback.isCanceled():
            return {}

        # Raster_θ - Incidence angle
        alg_params = {
            'EXTENT': parameters['dtm'],
            'NUMBER': parameters['__incidence_angle_090'],
            'OUTPUT_TYPE': 5,
            'PIXEL_SIZE': parameters['cell_size'],
            'TARGET_CRS': parameters['dtm'],
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['Raster_IncidenceAngle'] = processing.run('native:createconstantrasterlayer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(3)
        if feedback.isCanceled():
            return {}

        # H
        alg_params = {
            'BAND_A': 1,
            'BAND_B': None,
            'BAND_C': None,
            'BAND_D': None,
            'BAND_E': None,
            'BAND_F': None,
            'EXTENT_OPT': 0,
            'EXTRA': '',
            'FORMULA': 'cos(radians(A))',
            'INPUT_A': outputs['Raster_IncidenceAngle']['OUTPUT'],
            'INPUT_B': None,
            'INPUT_C': None,
            'INPUT_D': None,
            'INPUT_E': None,
            'INPUT_F': None,
            'NO_DATA': None,
            'OPTIONS': '',
            'PROJWIN': None,
            'RTYPE': 5,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['H'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(4)
        if feedback.isCanceled():
            return {}

        # N
        alg_params = {
            'BAND_A': 1,
            'BAND_B': 1,
            'BAND_C': None,
            'BAND_D': None,
            'BAND_E': None,
            'BAND_F': None,
            'EXTENT_OPT': 0,
            'EXTRA': '',
            'FORMULA': 'cos(radians(90)-radians(A))*cos(radians(180)-radians(B))',
            'INPUT_A': outputs['Raster_IncidenceAngle']['OUTPUT'],
            'INPUT_B': outputs['AzimutAngle']['OUTPUT'],
            'INPUT_C': None,
            'INPUT_D': None,
            'INPUT_E': None,
            'INPUT_F': None,
            'NO_DATA': None,
            'OPTIONS': '',
            'PROJWIN': None,
            'RTYPE': 5,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['N'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(5)
        if feedback.isCanceled():
            return {}

        #  S - Slope
        alg_params = {
            'AS_PERCENT': False,
            'BAND': 1,
            'COMPUTE_EDGES': False,
            'EXTRA': '',
            'INPUT': parameters['dtm'],
            'OPTIONS': '',
            'SCALE': 1,
            'ZEVENBERGEN': False,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['SSlope'] = processing.run('gdal:slope', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(6)
        if feedback.isCanceled():
            return {}

        # E
        alg_params = {
            'BAND_A': 1,
            'BAND_B': 1,
            'BAND_C': None,
            'BAND_D': None,
            'BAND_E': None,
            'BAND_F': None,
            'EXTENT_OPT': 0,
            'EXTRA': '',
            'FORMULA': 'cos(radians(90)-radians(A))*cos(radians(270)-radians(B))',
            'INPUT_A': outputs['Raster_IncidenceAngle']['OUTPUT'],
            'INPUT_B': outputs['AzimutAngle']['OUTPUT'],
            'INPUT_C': None,
            'INPUT_D': None,
            'INPUT_E': None,
            'INPUT_F': None,
            'NO_DATA': None,
            'OPTIONS': '',
            'PROJWIN': None,
            'RTYPE': 5,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['E'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(7)
        if feedback.isCanceled():
            return {}

        # A - Aspect
        alg_params = {
            'BAND': 1,
            'COMPUTE_EDGES': False,
            'EXTRA': '',
            'INPUT': parameters['dtm'],
            'OPTIONS': '',
            'TRIG_ANGLE': False,
            'ZERO_FLAT': False,
            'ZEVENBERGEN': False,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['AAspect'] = processing.run('gdal:aspect', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(8)
        if feedback.isCanceled():
            return {}

        # 1. Generazione C_Index
        alg_params = {
            'BAND_A': 1,
            'BAND_B': 1,
            'BAND_C': 1,
            'BAND_D': 1,
            'BAND_E': 1,
            'BAND_F': None,
            'EXTENT_OPT': 0,
            'EXTRA': '',
            'FORMULA': 'C * (cos(radians(A)))*(sin(radians(B)-radians(90))) + D * (-1) * (cos(radians(A)))*(cos((radians(B))-radians(90))) + E * (sin(radians(A)))',
            'INPUT_A': outputs['SSlope']['OUTPUT'],
            'INPUT_B': outputs['AAspect']['OUTPUT'],
            'INPUT_C': outputs['N']['OUTPUT'],
            'INPUT_D': outputs['E']['OUTPUT'],
            'INPUT_E': outputs['H']['OUTPUT'],
            'INPUT_F': None,
            'NO_DATA': None,
            'OPTIONS': '',
            'PROJWIN': parameters['dtm'],
            'RTYPE': 5,
            'OUTPUT': c_index_dest
        }
        outputs['C_index'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        
        # Forza il caricamento con il nome personalizzato
        context.addLayerToLoadOnCompletion(
            outputs['C_index']['OUTPUT'],
            QgsProcessingContext.LayerDetails(name_c_index, context.project(), 'C_index')
        )

        feedback.setCurrentStep(9)
        if feedback.isCanceled():
            return {}

        # 2. Generazione Percentage_c_index
        alg_params = {
            'BAND_A': 1,
            'BAND_B': None,
            'BAND_C': None,
            'BAND_D': None,
            'BAND_E': None,
            'BAND_F': None,
            'EXTENT_OPT': 0,
            'EXTRA': '',
            'FORMULA': ' abs(A)*100',
            'INPUT_A': outputs['C_index']['OUTPUT'],
            'INPUT_B': None,
            'INPUT_C': None,
            'INPUT_D': None,
            'INPUT_E': None,
            'INPUT_F': None,
            'NO_DATA': None,
            'OPTIONS': '',
            'PROJWIN': None,
            'RTYPE': 5,
            'OUTPUT': perc_c_index_dest
        }
        outputs['Percentage_c_index'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        
        # Forza il caricamento con il nome personalizzato
        context.addLayerToLoadOnCompletion(
            outputs['Percentage_c_index']['OUTPUT'],
            QgsProcessingContext.LayerDetails(name_percentage, context.project(), 'Percentage_c_index')
        )

        feedback.setCurrentStep(10)
        if feedback.isCanceled():
            return {}

        # 3. Generazione Class_of_visibility
        alg_params = {
            'DATA_TYPE': 5,
            'INPUT_RASTER': outputs['Percentage_c_index']['OUTPUT'],
            'NODATA_FOR_MISSING': False,
            'NO_DATA': -9999,
            'RANGE_BOUNDARIES': 0,
            'RASTER_BAND': 1,
            'TABLE': ['0','20','1','20','40','2','40','60','3','60','80','4','80','100','5'],
            'OUTPUT': class_vis_dest
        }
        outputs['Reclassify_percentage'] = processing.run('native:reclassifybytable', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        
        # Forza il caricamento con il nome personalizzato
        context.addLayerToLoadOnCompletion(
            outputs['Reclassify_percentage']['OUTPUT'],
            QgsProcessingContext.LayerDetails(name_visibility, context.project(), 'Class_of_visibility')
        )

        # --- APPLICAZIONE AUTOMATICA DELLO STILE TRAMITE POST-PROCESSOR ---
        plugin_dir = os.path.dirname(__file__)
        qml_path = os.path.join(plugin_dir, "stile_visibilita.qml")

        if context.willLoadLayerOnCompletion('Class_of_visibility'):
            global style_processor_instance
            style_processor_instance = VisibilityStyleProcessor(qml_path)
            context.layerToLoadOnCompletionDetails('Class_of_visibility').setPostProcessor(style_processor_instance)

        return {}

    def name(self):
        return 'Groundmotion - C index'

    def displayName(self):
        return 'Groundmotion - C index'

    def group(self):
        return 'Landslide'

    def groupId(self):
        return 'Landslide'

    def shortHelpString(self):
        return """<html><body><p><!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN" "http://www.w3.org/TR/REC-html40/strict.dtd">
<html><head><meta name="qrichtext" content="1" /><style type="text/css">
p, li { white-space: pre-wrap; }
</style></head><body style=" font-family:'MS Shell Dlg 2'; font-size:8.25pt; font-weight:400; font-style:normal;">
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-size:18pt;">C-INDEX</span></p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">The algorithm allows for the calculation of the C index, an indicator representing the percentage of detectable ground movement by a satellite. The C index ranges from -1 to 1 and is important for evaluating the capability of the satellite to monitor ground movements and identify deformations or terrain changes. In addition to the C-index map, the algorithm generates a map of the visibility percentage.This thematic map is useful for understanding how much of the terrain is actually visible from the satellite, taking into account the radar's acquisition angles. The visibility percentage helps determine the satellite's effectiveness in collecting accurate and comprehensive data. Finally, the algorithm produces a visibility class map. This map classifies different areas of the terrain into five classes based on their visibility, providing a detailed overview of the zones that can be monitored with greater precision compared to others.</p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">To balance accuracy and computational efficiency, we recommend a cell resolution equal to or greater than 20 m, according to InSAR resolution.</p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">The percentage of satellite-detectable movement can be represented by a map by calculating the coefficient C, which is derived from the following relationship: </p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-style:italic;">C = [N * cos(S) * sin (A - π/2)] + [E * (-1) * cos(S) * cos (A - π/2)] + [H * sin (S)]</span>   </p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">Where, </p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">S = Slope = slope of the side, calculated in radians and obtained from the digital soil model; </p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">A = Aspect = exposure of the side to the cardinal points, calculated</p></body></html>"""