Reservoir_algorithm.py

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

"""
/***************************************************************************
 Reservoir
                                 A QGIS plugin

 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2022-04-02
        copyright            : (C) 2022 by Iñigo Marin
        email                : proy.mbr@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__ = 'Iñigo Marin'
__date__ = '2022-04-02'
__copyright__ = '(C) 2022 by Iñigo Marin'

# 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
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterNumber
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsProcessingParameterPoint
from qgis.core import QgsProcessingParameterRasterDestination
from qgis.core import QgsProcessingParameterVectorDestination
from qgis.core import QgsProcessingParameterFeatureSink
from qgis.core import QgsExpression
from qgis.core import QgsMessageLog
from qgis.gui import QgsMessageBar
from qgis.core import Qgis
from os import path
from PyQt5.QtGui import QIcon
import processing


class ReservoirAlgorithm(QgsProcessingAlgorithm):
    """
    
    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.

    OUTPUT = 'OUTPUT'
    INPUT = 'INPUT'

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


        self.addParameter(QgsProcessingParameterRasterLayer('DEM', 'Existing DEM raster layer', defaultValue=None))
        self.addParameter(QgsProcessingParameterVectorLayer('dam_vector_layer', 'Dam vector layer', types=[QgsProcessing.TypeVector], defaultValue=None))
        self.addParameter(QgsProcessingParameterNumber('dam_top_level', 'Top dam level', type=QgsProcessingParameterNumber.Double, defaultValue=0))
        self.addParameter(QgsProcessingParameterNumber('water_depth_from_top', 'Water depth from top', type=QgsProcessingParameterNumber.Double, minValue=0, defaultValue=1))
        self.addParameter(QgsProcessingParameterPoint('reservoir_inner_point', 'Point inside the reservoir'))
        self.addParameter(QgsProcessingParameterNumber('dam_wall_slope', 'Dam wall slope 1:H', type=QgsProcessingParameterNumber.Double, minValue=0, defaultValue=3))
        self.addParameter(QgsProcessingParameterRasterDestination('new_lake_raster', 'New lake raster layer', createByDefault=True, defaultValue=''))
        self.addParameter(QgsProcessingParameterRasterDestination('new_DEM_raster', 'New DEM with dam raster layer', createByDefault=True, defaultValue=''))


    def processAlgorithm(self, parameters, context, model_feedback):
        """
        Here is where the processing itself takes place.
        """
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(6, model_feedback)
        results = {}
        outputs = {}

        # Búfer múlti-anillos (distancia constante)
        alg_params = {
            'DISTANCE': 1,
            'INPUT': parameters['dam_vector_layer'],
            'RINGS': 1000,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['BferMltianillosDistanciaConstante'] = processing.run('native:multiringconstantbuffer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

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

        # Calculadora de campos
        alg_params = {
            'FIELD_LENGTH': 10,
            'FIELD_NAME': 'top_level',
            'FIELD_PRECISION': 3,
            'FIELD_TYPE': 0,  # Float
            'FORMULA': parameters['dam_top_level'],
            'INPUT': outputs['BferMltianillosDistanciaConstante']['OUTPUT'],
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT 
        }
        outputs['CalculadoraDeCampos1'] = processing.run('native:fieldcalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
       
        # Calculadora de campos
        alg_params = {
            'FIELD_LENGTH': 10,
            'FIELD_NAME': 'wall_slope',
            'FIELD_PRECISION': 3,
            'FIELD_TYPE': 0,  # Float
            'FORMULA': parameters['dam_wall_slope'],
            'INPUT': outputs['CalculadoraDeCampos1']['OUTPUT'],
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['CalculadoraDeCampos2'] = processing.run('native:fieldcalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        
        # Advanced Python field calculator
        alg_params = {
            'FIELD_LENGTH': 10,
            'FIELD_NAME': 'fill_level',
            'FIELD_PRECISION': 3,
            'FIELD_TYPE': 1,  # Float
            'FORMULA': 'value = <top_level>- <distance>/<wall_slope>',
            'GLOBAL': '',
            'INPUT': outputs['CalculadoraDeCampos2']['OUTPUT'],
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['AdvancedPythonFieldCalculator'] = processing.run('qgis:advancedpythonfieldcalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

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

        # Rasterizar (vectorial a ráster)
        alg_params = {
            'BURN': None,
            'DATA_TYPE': 5,  # Float32
            'EXTENT': parameters['DEM'],
            'EXTRA': '',
            'FIELD': 'fill_level',
            'HEIGHT': 1,
            'INIT': 0,
            'INPUT': outputs['AdvancedPythonFieldCalculator']['OUTPUT'],
            'INVERT': False,
            'NODATA': None,
            'OPTIONS': '',
            'UNITS': 1,  # Georeferenced units
            'USE_Z': False,
            'WIDTH': 1,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        }
        outputs['RasterizarVectorialARster'] = processing.run('gdal:rasterize', alg_params, context=context, feedback=feedback, is_child_algorithm=True)


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

        # Raster calculator
        alg_params = {
            'BAND_A': 1,
            'BAND_B': 1,
            'BAND_C': None,
            'BAND_D': None,
            'BAND_E': None,
            'BAND_F': None,
            'EXTRA': '--hide',
            'FORMULA': '(A<B)*B + (A>=B)*A',
            'INPUT_A': parameters['DEM'],
            'INPUT_B': outputs['RasterizarVectorialARster']['OUTPUT'],
            'INPUT_C': None,
            'INPUT_D': None,
            'INPUT_E': None,
            'INPUT_F': None,
            'NO_DATA': None,
            'OPTIONS': '',
            'RTYPE': 5,  # Float32
            'OUTPUT': parameters['new_DEM_raster']
        }
        outputs['RasterCalculator'] = processing.run('gdal:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['new_DEM_raster'] = outputs['RasterCalculator']['OUTPUT']

        # Set layer style
        # Asigna el estilo de Pks
        
        alg_params = {
            'INPUT': outputs['RasterCalculator']['OUTPUT'],
            'STYLE': path.dirname(__file__) +'/DEM_Style.qml' 
        }              
        
        outputs['SetLayerStyle_DEM'] = processing.run('native:setlayerstyle', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        # r.lake
        alg_params = {
            '-n': False,
            'GRASS_RASTER_FORMAT_META': '',
            'GRASS_RASTER_FORMAT_OPT': '',
            'GRASS_REGION_CELLSIZE_PARAMETER': 0,
            'GRASS_REGION_PARAMETER': None,
            'coordinates': parameters['reservoir_inner_point'],
            'elevation': outputs['RasterCalculator']['OUTPUT'],
            'seed': None,
            'water_level': parameters['dam_top_level'] - parameters['water_depth_from_top'],
            'lake': parameters['new_lake_raster']
        }
        outputs['Rlake'] = processing.run('grass7:r.lake', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['new_lake_raster'] = outputs['Rlake']['lake']

        # Set layer style
        # Asigna el estilo
        
        alg_params = {
            'INPUT': outputs['Rlake']['lake'],
            'STYLE': path.dirname(__file__) +'/Lake_Style.qml' 
        }              
        
        outputs['SetLayerStyle_Lake'] = processing.run('native:setlayerstyle', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
   

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

        # Estadísticas de capa ráster
        alg_params = {
            'BAND': 1,
            'INPUT': outputs['Rlake']['lake']
        }
        outputs['EstadsticasDeCapaRster'] = processing.run('qgis:rasterlayerstatistics', alg_params, context=context, feedback=feedback, is_child_algorithm=True)


  
        return results

    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 'Reservoir'

    def icon(self):
        """
        Should return a QIcon which is used for your provider inside
        the Processing toolbox.
        """
        return QIcon(path.dirname(__file__) +'/Civileng.png')
    
    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 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 ''

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

    def shortHelpString(self):
        return """<html><body><p>Reservoir algorithm creates a dam on an existing DEM and calculates the flooded area</p>
<h2>Input parameters</h2>
<h3>Existing DEM raster layer</h3>
<p>Existing initial terrain DEM to be modified</p>
<h3>Dam vector layer</h3>
<p>Vector layer (line or polygon) containing just one line or polyline representing the crest of the dam</p>
<h3>Top dam level</h3>
<p>Absolute level of the dam crest</p>
<h3>Water depth from top</h3>
<p>Relative level of the water surface below the dam crest</p
<h3>Point inside the reservoir</h3>
<p>Point upstream of the dam wall inside the flooded area</p>
<h3>Dam wall slope</h3>
<p>Slope of the walls/enbankments of the dam as rate of vertical to horizontal lenght (1:H)</p>
<h2>Outputs</h2>
<h3>New lake raster layer</h3>
<p>New raster layer representing the water depth</p>
<h3>New DEM with dam raster layer</h3>
<p>Modified DEM raster layer with the dam</p>
<br><p align="right">Algorithm author: Iñigo Marin</p><p align="right">Algorithm version: v. 0.1</p></body></html>"""

    def helpUrl(self):
        return 'https://github.com/MBR111/CivilEng'
        
    def createInstance(self):
        return ReservoirAlgorithm()