Cindex_algorithm.py

# -*- 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$'

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterNumber
from qgis.core import QgsProcessingParameterRasterDestination
import processing


class GroundmotionCIndex(QgsProcessingAlgorithm):

    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):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(11, model_feedback)
        results = {}
        outputs = {}

        # 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,  # Ignore
            '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,  # Float32
            '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,  # Float32
            '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,  # Ignore
            '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,  # Float32
            '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,  # Ignore
            '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,  # Float32
            '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,  # Ignore
            '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,  # Float32
            '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 {}

        # 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,  # Ignore
            '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,  # Float32
            'OUTPUT': parameters['C_index']
        }
        outputs['C_index'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['C_index'] = outputs['C_index']['OUTPUT']

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

        # 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,  # Ignore
            '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,  # Float32
            'OUTPUT': parameters['Percentage_c_index']
        }
        outputs['Percentage_c_index'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Percentage_c_index'] = outputs['Percentage_c_index']['OUTPUT']

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

        # Reclassify_percentage
        alg_params = {
            'DATA_TYPE': 5,  # Float32
            'INPUT_RASTER': outputs['Percentage_c_index']['OUTPUT'],
            'NODATA_FOR_MISSING': False,
            'NO_DATA': -9999,
            'RANGE_BOUNDARIES': 0,  # min < valore <= max
            'RASTER_BAND': 1,
            'TABLE': ['0','20','1','20','40','2','40','60','3','60','80','4','80','100','5'],
            'OUTPUT': parameters['Class_of_visibility']
        }
        outputs['Reclassify_percentage'] = processing.run('native:reclassifybytable', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Class_of_visibility'] = outputs['Reclassify_percentage']['OUTPUT']
        return results

    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 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;">H, N, E = three cosines directors, the value of which depends on the satellite orbit parameters </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;">N = cos (π/2 - θ) * cos (π - φ)</span></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;">E = cos (π/2 - θ) * cos (3π/2 - φ)</span> </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;">H = cos (θ)</span> </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;">φ = track angle +90 = </span>azimuthal angle relative to satellite LOS </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;">θ =</span> incidence angle</p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">The track angle and incidence angle are provided within the processed SAR images for each satellite (and, in the case of EGMS, are found among the point attributes, while slope and aspect values are derived from the DTM.</p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">The aim is to calculate VSLOPE that is the projection of the LOS velocity (VLOS) along the direction of maximum slope through the relation:</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;">V</span><span style=" font-style:italic; vertical-align:sub;">slope</span><span style=" font-style:italic;"> = V</span><span style=" font-style:italic; vertical-align:sub;">los</span><span style=" font-style:italic;">/C</span> </p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">For C values close to 1, the movement is correctly observed by the satellite and the motion component is close to 100%. For negative <span style=" font-style:italic;">C</span> values, the recorded movement is in the opposite direction; in this case, if the speeds agree<span style=" color:#ff0000;"> </span>with the movement along the slope, the V<span style=" vertical-align:sub;">slope</span> will report a change of sign with respect to the original VLOS. For C values close to zero, the V<span style=" vertical-align:sub;">slope</span> tends to infinity, for this reason,<span style=" font-style:italic;"> C=0.2</span> is imposed when <span style=" font-style:italic;">-0.2&lt;C&lt;0.2</span>. </p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">In the analysis of landslide movements, the values recorded by the satellite may be significantly underestimated compared to the real values. Therefore, it is of fundamental importance to compare the data obtained from interferometric processing of satellite images with the percentage of movement detected along the LOS by that specific satellite.</p>
<p style=" margin-top:12px; margin-bottom:12px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;">Otherwise, one could wrongly attribute a &quot;stable&quot; state to a landslide movement without taking into account that the direction of observation of the satellite may not be suitable.</p>
<p style="-qt-paragraph-type:empty; margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><br /></p></body></html></p>
<h2>Parametri in ingresso
</h2>
<h3>DTM</h3>
<p>Raster layer of the DTM.
Preferable resolution greater than 20 m</p>
<h3>Cell size</h3>
<p>Dimension of the DTM cell (m)</p>
<h3>θ - Incidence angle (0-90)</h3>
<p>The inclination of the LOS with respect to the vertical.
Average value of 350° for the ascending orbit and 190° for the descending orbit.</p>
<h3>φ - Track angle - Azimut (0-360)</h3>
<p>The angle between the projection of the satellite trajectory on the horizontal plane and geographic north.
Average value of 39° for both orbits.</p>
<h2>Risultati</h2>
<h3>C_Index</h3>
<p>The result will be a raster in which the C-index values range from -1 to +1. This way you can identify the best or worst exposed slopes.</p>
<h3>Percentage_c_index</h3>
<p>Raster of the spatial distribution of the percentage of estimated real speed effectively detected along the ascending and descending LOS relative to the satellite SENTINEL-1.</p>
<h3>Class_of_visibility</h3>
<p>Classification of the percentage of visibility in 5 classes
Class 1: 0-20%
Class 2: 20-40%
Class 3: 40-60%
Class 4: 60-80%
Class 5: 80-100%</p>
<h2>Esempi</h2>
<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="-qt-paragraph-type:empty; margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><br /></p></body></html></p><br><p align="right">Autore della guida: Claudio Fasciano</p></body></html>"""

    def createInstance(self):
        return GroundmotionCIndex()