priorityWMCIFA.py

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

'''
/***************************************************************************
 DrainageBasinGeomorphology
                                 A QGIS plugin
 This plugin provides tools for geomorphological analysis in drainage basins.
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2025-03-22
        copyright            : (C) 2025 by João Vitor Pimenta
        email                : jvpjoaopimenta@gmail.com
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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__ = 'João Vitor Pimenta'
__date__ = '2025-03-22'
__copyright__ = '(C) 2025 by João Vitor Pimenta'

# This will get replaced with a git SHA1 when you do a git archive

__revision__ = '$Format:%H$'

import os
from qgis.PyQt.QtCore import QCoreApplication, QVariant
from qgis.PyQt.QtGui import QIcon
from qgis.core import (QgsProcessing,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterFileDestination,
                       QgsProcessingParameterNumber,
                       QgsProcessingParameterEnum,
                       QgsProcessingParameterBoolean,
                       QgsProcessingParameterFeatureSink,
                       QgsField)
from .algorithms.priorizationWMCIFA import calcWMCIFA,verifyLibs

class morphometricAnalysisWMCIFA(QgsProcessingAlgorithm):
    '''
    This is an example algorithm that takes a vector layer and
    creates a new identical one.

    It is meant to be used as an example of how to create your own
    algorithms and explain methods and variables used to do it. An
    algorithm like this will be available in all elements, and there
    is not need for additional work.

    All Processing algorithms should extend the QgsProcessingAlgorithm
    class.
    '''

    # Constants used to refer to parameters and outputs. They will be
    # used when calling the algorithm from another algorithm, or when
    # calling from the QGIS console.

    MORPHOMETRICS_PARAMETERS = 'MORPHOMETRICS_PARAMETERS'
    MORPHOMETRICS_PARAMETERS_STANDARD = 'MORPHOMETRICS_PARAMETERS_STANDARDIZED'
    CORRELATION_MATRIX = 'CORRELATION_MATRIX'
    TOTAL_VARIANCE_EXPLAINED_TABLE = 'TOTAL_VARIANCE_EXPLAINED_TABLE'
    UNROTATED_AND_ROTATED_MATRIX = 'UNROTATED_AND_ROTATED_MATRIX'
    RANKING_TABLE_WITH_CP_VALUES = 'RANKING_TABLE_WITH_CP_VALUES'
    DRAINAGE_BASINS = 'DRAINAGE_BASINS'
    DEM = 'DEM'
    CHANNEL_COORDINATE_PRECISION = 'CHANNEL_COORDINATE_PRECISION'
    CHANNEL_NETWORK = 'CHANNEL_NETWORK'
    SELECTED_PARAMETERS_DIRECTLY_PROPORTIONAL = 'SELECTED_PARAMETERS_DIRECTLY_PROPORTIONAL'
    SELECTED_PARAMETERS_INVERSELY_PROPORTIONAL = 'SELECTED_PARAMETERS_INVERSELY_PROPORTIONAL'
    BASINS_RANKED = 'BASINS_RANKED'
    DECIMAL_PLACES = 'DECIMAL_PLACES'
    MINIMUM_CHANNEL_LENGTH = 'MINIMUM_CHANNEL_LENGTH'
    LIMIT_FOR_VALLEY_FLOOR = 'LIMIT_FOR_VALLEY_FLOOR'
    MIN_FOR_VALLEY_HEIGHT = 'MIN_FOR_VALLEY_HEIGHT'
    USE_LONGEST_DRAINAGE = 'USE_LONGEST_DRAINAGE'
    POINTS_TTSF = 'POINTS_TTSF'
    N_SECTIONS_SL = 'N_SECTIONS_SL'
    POINTS_MIDLINE = 'POINTS_MIDLINE'


    def initAlgorithm(self, config):
        '''
        Here we define the inputs and output of the algorithm, along
        with some other properties.
        '''

        # We add the input vector features source. It can have any kind of
        # geometry.
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.DRAINAGE_BASINS,
                self.tr('Drainage basins'),
                [QgsProcessing.SourceType.TypeVectorPolygon]
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.CHANNEL_NETWORK,
                self.tr('Channel network'),
                [QgsProcessing.SourceType.TypeVectorLine]
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.DEM,
                self.tr('DEM'),
                [QgsProcessing.SourceType.TypeRaster]
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.N_SECTIONS_SL,
                self.tr('Number of sections for SL'),
                type=QgsProcessingParameterNumber.Type.Integer,
                minValue=1,
                defaultValue=10
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.LIMIT_FOR_VALLEY_FLOOR,
                self.tr('Limit for valley floor'),
                type=QgsProcessingParameterNumber.Type.Double,
                minValue=0,
                defaultValue=0.1,
                optional=False
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.MIN_FOR_VALLEY_HEIGHT,
                self.tr('Minimum height for valley peak'),
                type=QgsProcessingParameterNumber.Type.Double,
                minValue=0,
                defaultValue=1.0,
                optional=False
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.POINTS_MIDLINE,
                self.tr('Number of points to create midline'),
                type=QgsProcessingParameterNumber.Type.Integer,
                minValue=2,
                defaultValue=50,
                optional=False
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.POINTS_TTSF,
                self.tr('Number of points for TTSF'),
                type=QgsProcessingParameterNumber.Type.Integer,
                minValue=0,
                defaultValue=50,
                optional=False
            )
        )

        self.addParameter(
            QgsProcessingParameterBoolean(
                self.USE_LONGEST_DRAINAGE,
                self.tr('Use lch as longest drainage and not the main channel'),
                defaultValue=False
            )
        )

        self.parametersToChoose = ['Mean stream length ratio',
                         'Mean bifurcation ratio',
                         'RHO coefficient',
                         'Main channel sinuosity index',
                         'Fitness ratio (Rf)',
                         'Wandering ratio (Rw)',
                         'Drainage density (Dd) (km/km2)',
                         'Stream frequency (Fs) (1/km2)',
                         'Drainage texture (Dt) (1/km)',
                         'Length of overland flow (Lo) (km)',
                         'Constant of channel maintenance (Ccm) (km2/km)',
                         'Drainage intensity (Di) (1/km)',
                         'Infiltration number (If) (km/km4)',
                         'Area (km2)',
                         'Perimeter (km)',
                         'Basin length (Lg) (km)',
                         'Circulatory ratio (Rc)',
                         'Elongation ratio (Re)',
                         'Form factor (Ff)',
                         'Lemniscate ratio (K)',
                         'Shape index (Sb)',
                         'Compactness coefficient (Cc)',
                         'Minimum elevation (m)',
                         'Maximum elevation (m)',
                         'Mean elevation (m)',
                         'Relief (Bh) (m)',
                         'Relief ratio (Rh)',
                         'Relative relief (Rhp)',
                         'Ruggedness number (Rn)',
                         'Dissection index (Di)',
                         'Gradient ratio (Gr)',
                         'Transverse topographic symmetry factor (TTSF)',
                         'Assimetry factor (AF)',
                         'Stream-Length index mean (SLm)',
                         'Stream-Length index total mean (SLtm)',
                         'Stream-Length index (SLm/SLtm)',
                         'Valley floor width-height ratio (Vf)',
                         'None']

        self.addParameter(
            QgsProcessingParameterEnum(
                self.SELECTED_PARAMETERS_DIRECTLY_PROPORTIONAL,
                description='Parameters for WMCI-FA analysis (directly proportional)',
                options=self.parametersToChoose,
                allowMultiple=True,
                defaultValue=[],
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.SELECTED_PARAMETERS_INVERSELY_PROPORTIONAL,
                description='Parameters for WMCI-FA analysis (inversely proportional)',
                options=self.parametersToChoose,
                allowMultiple=True,
                defaultValue=[],
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.CHANNEL_COORDINATE_PRECISION,
                self.tr('Channel coordinate precision to snap'),
                type=QgsProcessingParameterNumber.Type.Double,
                minValue=0,
                defaultValue=0.01,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.MINIMUM_CHANNEL_LENGTH,
                self.tr('Minimum channel length'),
                type=QgsProcessingParameterNumber.Type.Double,
                minValue=0,
                defaultValue=0.01,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.DECIMAL_PLACES,
                self.tr('Decimal places of the result'),
                type=QgsProcessingParameterNumber.Type.Integer,
                minValue=0,
                defaultValue=2,
                optional=False
            )
        )

        # We add a feature sink in which to store our processed features (this
        # usually takes the form of a newly created vector layer when the
        # algorithm is run in QGIS).
        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.MORPHOMETRICS_PARAMETERS,
                self.tr('Morphometric parameters'),
                fileFilter=('CSV files (*.csv)')
            )
        )

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.MORPHOMETRICS_PARAMETERS_STANDARD,
                self.tr('Morphometric parameters standardized'),
                fileFilter=('CSV files (*.csv)')
            )
        )

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.CORRELATION_MATRIX,
                self.tr('Correlation table'),
                fileFilter=('CSV files (*.csv)')
            )
        )

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.TOTAL_VARIANCE_EXPLAINED_TABLE,
                self.tr('Total variance explained table'),
                fileFilter=('CSV files (*.csv)')
            )
        )

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.UNROTATED_AND_ROTATED_MATRIX,
                self.tr('Unrotated and rotated matrix'),
                fileFilter=('CSV files (*.csv)')
            )
        )

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.RANKING_TABLE_WITH_CP_VALUES,
                self.tr('Ranking table with compound parameter values'),
                fileFilter=('CSV files (*.csv)')
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.BASINS_RANKED,
                self.tr('Ranked basins WMCI-FA'))
            )

    def processAlgorithm(self, parameters, context, feedback):
        '''
        Here is where the processing itself takes place.
        '''

        # Retrieve the feature source and sink. The 'dest_id' variable is used
        # to uniquely identify the feature sink, and must be included in the
        # dictionary returned by the processAlgorithm function.
        basinSource = self.parameterAsSource(parameters, self.DRAINAGE_BASINS, context)

        channelNetwork = self.parameterAsSource(parameters, self.CHANNEL_NETWORK, context)

        demLayer = self.parameterAsRasterLayer(parameters, self.DEM, context)

        pointsTTSF = self.parameterAsInt(parameters, self.POINTS_TTSF, context)

        limitForValleyFloor = self.parameterAsDouble(parameters, self.LIMIT_FOR_VALLEY_FLOOR, context)

        minForValleyHeight = self.parameterAsDouble(parameters, self.MIN_FOR_VALLEY_HEIGHT, context)

        useLongestDrainage = self.parameterAsBoolean(parameters, self.USE_LONGEST_DRAINAGE, context)

        precisionSnapCoordinates = self.parameterAsDouble(parameters, self.CHANNEL_COORDINATE_PRECISION, context)
        minimumChannelLength = self.parameterAsDouble(parameters, self.MINIMUM_CHANNEL_LENGTH, context)

        selectedParametersDirectly = self.parameterAsEnums(parameters, self.SELECTED_PARAMETERS_DIRECTLY_PROPORTIONAL, context)
        selectedStringsDirectly = [self.parametersToChoose[i] for i in selectedParametersDirectly]

        selectedParametersInversely = self.parameterAsEnums(parameters, self.SELECTED_PARAMETERS_INVERSELY_PROPORTIONAL, context)
        selectedStringsInversely = [self.parametersToChoose[i] for i in selectedParametersInversely]

        decimalPlaces = self.parameterAsInt(parameters, self.DECIMAL_PLACES, context)

        pathParameters = self.parameterAsFileOutput(parameters, self.MORPHOMETRICS_PARAMETERS, context)
        pathParametersStandardized = self.parameterAsFileOutput(parameters, self.MORPHOMETRICS_PARAMETERS_STANDARD, context)
        pathCorrMatrix = self.parameterAsFileOutput(parameters, self.CORRELATION_MATRIX, context)
        pathVarExplained = self.parameterAsFileOutput(parameters, self.TOTAL_VARIANCE_EXPLAINED_TABLE, context)
        pathRotUnrot = self.parameterAsFileOutput(parameters, self.UNROTATED_AND_ROTATED_MATRIX, context)
        pathRankCp = self.parameterAsFileOutput(parameters, self.RANKING_TABLE_WITH_CP_VALUES, context)

        nSectionsSL = self.parameterAsInt(parameters, self.N_SECTIONS_SL, context)
        pointsMidline = self.parameterAsInt(parameters, self.POINTS_MIDLINE, context)

        fields = basinSource.fields()
        fields.append(QgsField("ranking", QVariant.Double))
        fields.append(QgsField("priority", QVariant.String))

        basinsRanked, destId = self.parameterAsSink(
            parameters,
            self.BASINS_RANKED,
            context,
            fields,
            basinSource.wkbType(),
            basinSource.sourceCrs()
        )

        verifyLibs()
        calcWMCIFA(basinSource,channelNetwork,demLayer,feedback,precisionSnapCoordinates,decimalPlaces,selectedStringsDirectly,selectedStringsInversely,pathCorrMatrix,pathVarExplained,pathRotUnrot,pathRankCp,basinsRanked,pathParameters,minimumChannelLength,pathParametersStandardized,pointsTTSF,limitForValleyFloor,minForValleyHeight,useLongestDrainage,nSectionsSL,pointsMidline)

        # Return the results of the algorithm. In this case our only result is
        # the feature sink which contains the processed features, but some
        # algorithms may return multiple feature sinks, calculated numeric
        # statistics, etc. These should all be included in the returned
        # dictionary, with keys matching the feature corresponding parameter
        # or output names.
        return {self.MORPHOMETRICS_PARAMETERS: pathParameters,
                self.MORPHOMETRICS_PARAMETERS_STANDARD: pathParametersStandardized,
                self.CORRELATION_MATRIX: pathCorrMatrix,
                self.TOTAL_VARIANCE_EXPLAINED_TABLE: pathVarExplained,
                self.UNROTATED_AND_ROTATED_MATRIX: pathRotUnrot,
                self.RANKING_TABLE_WITH_CP_VALUES: pathRankCp,
                self.BASINS_RANKED: destId}

    def name(self):
        '''
        Returns the algorithm name, used for identifying the algorithm. This
        string should be fixed for the algorithm, and must not be localised.
        The name should be unique within each provider. Names should contain
        lowercase alphanumeric characters only and no spaces or other
        formatting characters.
        '''
        return 'Calculate basin priority (WMCI-FA)'

    def displayName(self):
        '''
        Returns the translated algorithm name, which should be used for any
        user-visible display of the algorithm name.
        '''
        return self.tr(self.name())

    def groupId(self):
        return "basin_priority"

    def group(self):
        '''
        Returns the name of the group this algorithm belongs to. This string
        should be localised.
        '''
        return self.tr("Basin priority")

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

    def shortHelpString(self):
        """
        Returns a localised short help string for the algorithm.
        """
        return self.tr("""
        <html>
            <body>
                <p>       
        This tool calculates all morphometric parameters of each basin feature individually and then calculates a priority order for the basins, based on the parameters selected by the user and the morphometric method.               
                </p>
                <p>
        <strong>Drainage basins: </strong>Layer containing drainage basins as features.
        <strong>Channel network: </strong>Layer containing the drainage network of the drainage basins.
        <strong>DEM: </strong>Raster containing the band with the altimetry of the drainage basins.
        <strong>Number of sections for SL index: </strong>It is the number of sections used to calculate the SL index.
        <strong>Limit for valley floor: </strong>Its the height limit to calculate the valley floor.
        <strong>Minimum height for valley: </strong>Its the minimum height difference between one point and and next point to consider a valley peak. A valley peak is defined as the first point preceding a downward slope. A minimum threshold is applied to prevent minor dips which could result from inaccuracies from being counted.
        <strong>Number of points for midline: </strong>Its the number of points to create the midline of the basin.
        <strong>Number of points for TTSF: </strong>Its the number of points in the midline to calculate the TTSF.
        <strong>Use lch as longest drainage and not the main channel: </strong>The plugin default is to use the lch as main channel (hightest strahler order) but in some cases lch as longest drainage can be more useful. If this box is checked, the largest channel will be used as LCH in the calculations (see README).
        <strong>Parameters for WMCI-FA analysis (directly proportional): </strong>Morphometric parameters directly proportional to the priority the user wants to analyze.
        <strong>Parameters for WMCI-FA analysis (indirectly proportional): </strong>Morphometric parameters indirectly proportional to the priority the user wants to analyze.
        <strong>Channel coordinate precision: </strong>It is the precision of the channel coordinates, for example: for a precision of 0.000001 the coordinate xxxxxx.xxxxxxxxxxxx becomes xxxxxx.xxxxxx. It is recommended to use 0.000001 to correct possible geometry errors when selecting channels that intersect the basin. If it is 0, there will be no rounding.
        <strong>Minimum channel length: </strong>It is used to correct intersection errors, as well as channel network precision.
        <strong>Decimal places of the result: </strong>Number of decimal places in results.
        <strong>Morphometric parameters: </strong>File with all morphometric parameters calculated individually for each basin.
        <strong>Correlation table: </strong>File containing the table with the intercorrelation matrix of the parameters selected by the user.
        <strong>Total variance explained table: </strong>Amount of variance explained by each component of the PCA analysis (eigenvectors and eigenvalues), selected based on Kaiser's rule (eigenvectors with eigenvalues > 1).
        <strong>Unrotated and rotated matrix: </strong>Matrix with the correlation of the parameters in relation to the selected components, unrotated and rotated by the varimax method.
        <strong>Ranking table with compound parameter values: </strong>Table with ranked, weighted parameters and compound values with the final ranking.
        <strong>Ranked basins WMCI-FA: </strong>Original basin vector with two columns added, one with the compound parameter and the other with the final basin ranking.
                       
        The use of a projected CRS is recommended (the plugin calculation assumes that all input layers are in projected coordinate reference systems).
                       
        If you need more information about how the plugin works, such as the calculations it performs, among other things, access: https://github.com/JoaoVitorPimenta/qgis-plugin-Drainage-Basin-Geomorphology
        If you have found any bugs, errors or have any requests to make, among other things, please acess: https://github.com/JoaoVitorPimenta/qgis-plugin-Drainage-Basin-Geomorphology/issues
        If you need training for the plugin, or want to contact the plugin author for any reason, send an email to: jvpjoaopimentadev@gmail.com                </p>
            </body>
        </html>
                    """)

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return morphometricAnalysisWMCIFA()