runCom9MoTVoellmy_algorithm.py

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

"""
/***************************************************************************
 AvaFrameRunCom1DFA
                                 A QGIS plugin
 Connects to AvaFrame
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2021-08-26
        copyright            : (C) 2021 by AvaFrame Team
        email                : felix@avaframe.org
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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__ = "AvaFrame Team"
__date__ = "2022"
__copyright__ = "(C) 2022 by AvaFrame Team"

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

__revision__ = "$Format:%H$"


from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (
    QgsProcessing,
    QgsProcessingException,
    QgsProcessingAlgorithm,
    QgsProcessingParameterRasterLayer,
    QgsProcessingParameterMultipleLayers,
    QgsProcessingParameterFolderDestination,
    QgsProcessingParameterEnum,
    QgsProcessingOutputVectorLayer,
    QgsProcessingOutputMultipleLayers,
)


class runCom9MoTVoellmyAlgorithm(QgsProcessingAlgorithm):
    """
    This is the AvaFrame Connection, i.e. the part running with QGis. For this
    connector to work, more installation is needed. See instructions at docs.avaframe.org
    """

    DEM = "DEM"
    RELSHP = "RELSHP"
    RELRAS = "RELRAS"
    FRICTION = "FRICTION"
    MU = "MU"
    K = "K"
    ENTRAINMENT = "ENTRAINMENT"
    B0 = "B0"
    TAUC = "TAUC"
    FOREST = "FOREST"
    ND = "ND"
    BHD = "BHD"
    OUTPUT = "OUTPUT"
    OUTPPR = "OUTPPR"
    FOLDEST = "FOLDEST"

    ADDTONAME = "ADDTONAME"
    SMALLAVA = "SMALLAVA"
    DATA_TYPE = "DATA_TYPE"

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

        self.addParameter(
            QgsProcessingParameterRasterLayer(self.DEM, self.tr("DEM layer (.asc raster)"))
        )

        self.addParameter(
            QgsProcessingParameterMultipleLayers(
                self.RELSHP,
                self.tr(
                    "EITHER Release layer(s) as shape (with thickness attribute for release height)"
                ),
                layerType=QgsProcessing.TypeVectorAnyGeometry,
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterMultipleLayers(
                self.RELRAS,
                self.tr("OR Release layer(s) as .asc raster"),
                layerType=QgsProcessing.TypeRaster,
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.FRICTION,
                self.tr("-----------Friction-----------"),
                options=[
                    self.tr("Constant"),
                    self.tr("Variable (needs mu and k)"),
                ],
                defaultValue=0,
                allowMultiple=False,
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.MU,
                self.tr("mu / Dry friction coeff (.asc raster)"),
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.K,
                self.tr("k / Turbulent friction coeff (.asc raster)"),
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.ENTRAINMENT,
                self.tr("---------Entrainment----------"),
                options=[
                    self.tr("No Entrainment"),
                    self.tr("TJEM (needs b0 and tau_c)"),
                ],
                defaultValue=0,
                allowMultiple=False,
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.B0,
                self.tr("b0 / Erodible snow depth (.asc raster)"),
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.TAUC,
                self.tr("tau_c / Snow shear strength (.asc raster)"),
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.FOREST,
                self.tr("------------Forest------------"),
                options=[
                    self.tr("No Forest"),
                    self.tr("Forest (needs nd and bhd)"),
                ],
                defaultValue=0,
                allowMultiple=False,
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.ND,
                self.tr("nd / Forest density (.asc raster)"),
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.BHD,
                self.tr("bhd / Tree diameter (.asc raster)"),
                optional=True,
                defaultValue="",
            )
        )

        self.addParameter(
            QgsProcessingParameterFolderDestination(
                self.FOLDEST, self.tr("------Destination folder (Should be empty)------")
            )
        )

        self.addOutput(
            QgsProcessingOutputVectorLayer(
                self.OUTPUT,
                self.tr("Output layer"),
                QgsProcessing.TypeVectorAnyGeometry,
            )
        )

        self.addOutput(
            QgsProcessingOutputMultipleLayers(
                self.OUTPPR,
            )
        )

    def flags(self):
        return super().flags()
        # return super().flags() | QgsProcessingAlgorithm.FlagNoThreading

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

        import avaframe.version as gv
        from . import avaframeConnector_commonFunc as cF

        feedback.pushInfo("AvaFrame Version: " + gv.getVersion())

        targetADDTONAME = ""

        sourceDEM = self.parameterAsRasterLayer(parameters, self.DEM, context)
        if sourceDEM is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.DEM))

        # Release files - check both shapefile and raster options
        allRELSHP = self.parameterAsLayerList(parameters, self.RELSHP, context)
        allRELRAS = self.parameterAsLayerList(parameters, self.RELRAS, context)

        # Check that only one release type is provided
        if allRELSHP and allRELRAS:
            raise QgsProcessingException(
                self.tr(
                    "Error: Please provide EITHER release layers as shapefile OR as raster, not both"
                )
            )

        relShpDict = {}
        if allRELSHP:
            relShpDict = {lyr.source(): lyr for lyr in allRELSHP}

        relRasDict = {}
        if allRELRAS:
            relRasDict = {lyr.source(): lyr for lyr in allRELRAS}

        # Get friction option
        frictionOption = self.parameterAsInt(parameters, self.FRICTION, context)
        frictionOptions = ["constant", "variable"]
        frictionString = frictionOptions[frictionOption]

        # Extract raster parameters
        sourceMU = self.parameterAsRasterLayer(parameters, self.MU, context)
        sourceK = self.parameterAsRasterLayer(parameters, self.K, context)

        # Validate friction-related parameters
        if frictionString == "variable":
            # Variable friction requires MU and K
            if sourceMU is None:
                raise QgsProcessingException(
                    self.tr("Error: Variable friction requires mu (dry friction coeff) to be provided")
                )
            if sourceK is None:
                raise QgsProcessingException(
                    self.tr("Error: Variable friction requires k (turbulent friction coeff) to be provided")
                )

        # Get entrainment option
        entrainmentOption = self.parameterAsInt(parameters, self.ENTRAINMENT, context)
        entrainmentOptions = ["noEntrainment", "tjem"]
        entrainmentString = entrainmentOptions[entrainmentOption]

        sourceB0 = self.parameterAsRasterLayer(parameters, self.B0, context)
        sourceTAUC = self.parameterAsRasterLayer(parameters, self.TAUC, context)

        # Validate entrainment-related parameters
        hasEntrainment = 0
        if entrainmentString == "tjem":
            # TJEM requires K and B0
            if sourceTAUC is None:
                raise QgsProcessingException(
                    self.tr("Error: TJEM entrainment requires tau_c (snow shear strength ) to be provided")
                )
            if sourceB0 is None:
                raise QgsProcessingException(
                    self.tr("Error: TJEM entrainment requires b0 (erodible snow depth) to be provided")
                )
            hasEntrainment = 1

        # Get forest option
        forestOption = self.parameterAsInt(parameters, self.FOREST, context)
        forestOptions = ["noForest", "forest"]
        forestString = forestOptions[forestOption]

        sourceND = self.parameterAsRasterLayer(parameters, self.ND, context)
        sourceBHD = self.parameterAsRasterLayer(parameters, self.BHD, context)

        # Validate forest-related parameters
        hasForest = 0
        if forestString == "forest":
            # Forest requires ND and BHD
            if sourceND is None:
                raise QgsProcessingException(
                    self.tr("Error: Forest requires nd (forest density) to be provided")
                )
            if sourceBHD is None:
                raise QgsProcessingException(
                    self.tr("Error: Forest requires bhd (tree diameter) to be provided")
                )
            hasForest = 1

        sourceFOLDEST = self.parameterAsFile(parameters, self.FOLDEST, context)

        # create folder structure (targetDir is the tmp one)
        finalTargetDir, targetDir = cF.createFolderStructure(sourceFOLDEST)

        feedback.pushInfo(sourceDEM.source())

        # copy DEM
        cF.copyDEM(sourceDEM, targetDir)

        # copy all release shapefile parts
        if relShpDict:
            cF.copyMultipleShp(
                relShpDict, targetDir / "Inputs" / "REL", targetADDTONAME
            )

        # copy all release raster parts
        if relRasDict:
            for source in relRasDict:
                cF.copyRaster(relRasDict[source], targetDir / "Inputs" / "REL", "")

        # copy raster parameter files to RASTERS folder with suffixes
        if sourceMU is not None:
            cF.copyRaster(sourceMU, targetDir / "Inputs" / "RASTERS", "_mu")

        if sourceK is not None:
            cF.copyRaster(sourceK, targetDir / "Inputs" / "RASTERS", "_k")

        if sourceB0 is not None:
            cF.copyRaster(sourceB0, targetDir / "Inputs" / "ENT", "")

        if sourceTAUC is not None:
            cF.copyRaster(sourceTAUC, targetDir / "Inputs" / "RASTERS", "_tauc")

        if sourceND is not None:
            cF.copyRaster(sourceND, targetDir / "Inputs" / "RES", "_nd")

        if sourceBHD is not None:
            cF.copyRaster(sourceBHD, targetDir / "Inputs" / "RASTERS", "_bhd")

        feedback.pushInfo("Starting the simulations")
        feedback.pushInfo("This might take a while")
        feedback.pushInfo("See console for progress")

        # Determine -st argument based on FOREST and ENTRAINMENT parameters
        if hasForest and hasEntrainment:
            stValue = "entres"
        elif hasForest:
            stValue = "res"
        elif hasEntrainment:
            stValue = "ent"
        else:
            stValue = "null"

        # Generate command and run via subprocess.run
        command = [
            "python",
            "-m",
            "avaframe.com9MoTVoellmy.runCom9MoTVoellmy",
            str(targetDir),
            "-st",
            stValue,
        ]
        cF.runAndCheck(command, self, feedback)

        feedback.pushInfo("Done, start loading the results")

        # Move input, log and output folders to finalTargetDir
        cF.moveInputAndOutputFoldersToFinal(targetDir, finalTargetDir)

        # Get peakfiles to return to QGIS
        try:
            rasterResults = cF.getLatestPeakCom9(finalTargetDir)
        except:
            raise QgsProcessingException(
                self.tr(
                    "Something went wrong with com9MoTVoellmy, please check log files"
                )
            )

        allRasterLayers = cF.addStyleToCom1DFAResults(rasterResults)

        context = cF.addLayersToContext(context, allRasterLayers, self.OUTPPR)

        feedback.pushInfo("\n---------------------------------")
        feedback.pushInfo("Done, find results and logs here:")
        feedback.pushInfo(str(finalTargetDir.resolve()))
        feedback.pushInfo("---------------------------------\n")

        return {self.OUTPPR: allRasterLayers}

    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 "com9motvoellmy"

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

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

    def groupId(self):
        """
        Returns the unique ID of the group this algorithm belongs to. This
        string should be fixed for the algorithm, and must not be localised.
        The group id should be unique within each provider. Group id should
        contain lowercase alphanumeric characters only and no spaces or other
        formatting characters.
        """
        return "NGI_Experimental"

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

    def shortHelpString(self) -> str:
        hstring = "Runs dense flow simulations via module com9MoTVoellmy. \n\
                All raster files need to have the same extend as the DEM. \n\
                For more information go to (or use the help button below): \n\
                AvaFrame Documentation: https://docs.avaframe.org\n\
                Homepage: https://avaframe.org\n\
                Praxisleitfaden: https://avaframe.org/reports\n"

        return self.tr(hstring)
        # Praxisleitfaden: https://info.bml.gv.at/dam/jcr:edebd872-2a86-4edf-ac5e-635ef11e35fe/Praxisleitfaden%20LawSim%20WLV%202022%20Gr%C3%BCn.pdf\n'

    def helpUrl(self):
        return "https://docs.avaframe.org/en/latest/connector.html"

    def createInstance(self):
        return runCom9MoTVoellmyAlgorithm()