map_the_value_of_individual_ecosystem_services.py

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

"""
/***************************************************************************
 EcoValuator
                                 A QGIS plugin
 Calculate ecosystem service values for a given area
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2018-04-02
        copyright            : (C) 2018 by Key-Log Economics
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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__ = 'Key-Log Economics'
__date__ = '2018-04-02'
__copyright__ = '(C) 2018 by Key-Log Economics'

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

__revision__ = '$Format:%H$'

import processing

from os.path import splitext

from PyQt5.QtGui import *
from qgis.utils import *
from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterRasterDestination,
                       QgsRasterFileWriter,
                       QgsRasterLayer,
                       QgsProcessingParameterString,
                       QgsProcessingParameterNumber,
                       QgsProcessingOutputLayerDefinition,
                       QgsProcessingParameterEnum,
                       QgsMapLayerStyle,
                       QgsProject,
                       QgsColorRampShader,
                       QgsRasterShader,
                       QgsRasterBandStats,
                       QgsSingleBandPseudoColorRenderer,
                       QgsClassificationQuantile
                      )

import numpy as np

from .appinter import (Raster, App)
from .eco_valuator_classes import LULC_dataset, ESV_dataset, Symbology



class MapTheValueOfIndividualEcosystemServices(QgsProcessingAlgorithm):
    """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."""
    INPUT_RASTER = 'INPUT_RASTER'
    INPUT_ESV_FIELD = 'INPUT_ESV_FIELD'
    INPUT_ESV_STAT = 'INPUT_ESV_STAT'
    STATS_MAP = {'Minimum':'Min', 'Average':'Avg', 'Maximum':'Max'}
    STATS = list(STATS_MAP)
    OUTPUT_RASTER = 'OUTPUT_RASTER'
    OUTPUT_RASTER_FILENAME_DEFAULT = 'Output esv raster'
    INPUT_LULC_SOURCE = 'INPUT_LULC_SOURCE'

    with ESV_dataset() as esv:
        LULC_SOURCES = esv.get_lulc_sources()
        INPUT_ESV_FIELD_OPTIONS = esv.get_ecosystem_service_names()

    def initAlgorithm(self, config):
        """
        Here we define the inputs and output of the algorithm
        """
        # Add a parameter for the clipped raster layer
        self.addParameter(
            QgsProcessingParameterEnum(
                self.INPUT_LULC_SOURCE,
                self.tr('Select land use/land cover data source'),
                self.LULC_SOURCES
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INPUT_RASTER,
                self.tr('Input clipped raster layer')
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.INPUT_ESV_FIELD,
                self.tr('Choose ecosystem service of interest'),
                self.INPUT_ESV_FIELD_OPTIONS
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.INPUT_ESV_STAT,
                self.tr('Choose ecosystem Service Value Level'),
                self.STATS
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterDestination(
                self.OUTPUT_RASTER,
                self.tr(self.OUTPUT_RASTER_FILENAME_DEFAULT)
            )
        )

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

        log = feedback.setProgressText
        input_lulc_source_index = self.parameterAsEnum(parameters, self.INPUT_LULC_SOURCE, context)
        input_lulc_source = self.LULC_SOURCES[input_lulc_source_index]
        input_raster = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER, context)
        input_esv_field_index = self.parameterAsEnum(parameters, self.INPUT_ESV_FIELD, context)
        input_esv_field = self.INPUT_ESV_FIELD_OPTIONS[input_esv_field_index]
        input_esv_stat_index = self.parameterAsEnum(parameters, self.INPUT_ESV_STAT, context)
        input_esv_stat_full_name = self.STATS[input_esv_stat_index]
        input_esv_stat = self.STATS_MAP[input_esv_stat_full_name]
        output_raster_destination = self.parameterAsOutputLayer(parameters, self.OUTPUT_RASTER, context)
        result = {self.OUTPUT_RASTER: output_raster_destination}


        #  // STEP 1. Check output file format to make sure it is a geotiff //
        output_format = QgsRasterFileWriter.driverForExtension(splitext(output_raster_destination)[1])

        if not output_format or output_format.lower() != "gtiff":
            error_message = "CRITICAL: Currently only GeoTIFF output format allowed, exiting!"
            feedback.reportError(error_message)
            return({'error': error_message})
        else:
            message = "Output file is GeoTIFF. Check"
            log(message)




        #  // STEP 2. Make instance of LULC dataset from clipped layer here //
        LULC_raster = LULC_dataset(input_lulc_source, input_raster)
        
        # Check to make sure all land use codes are valid 
        valid = LULC_raster.is_valid()
        if isinstance(valid, str):
            #If is instance returns a string it is not valid. The string contains the error message
            error_message = valid
            feedback.reportError(error_message)
            return {'error': error_message}




        #  // STEP 3. Reclassification //
        # Get reclassify table for selected parameters
        ESV_data = ESV_dataset()
        reclass_table = ESV_data.make_reclassify_table(LULC_raster.cell_size(),
                                                       input_lulc_source,
                                                       input_esv_stat,
                                                       input_esv_field)

        # Perform reclassification
        reclassify_params = {'INPUT_RASTER':input_raster,
                             'RASTER_BAND':1,
                             'TABLE':reclass_table,
                             'NO_DATA':-9999,
                             'RANGE_BOUNDARIES':0,
                             'NODATA_FOR_MISSING':True,
                             'DATA_TYPE':6,
                             'OUTPUT':output_raster_destination}

        processing.run("native:reclassifybytable", reclassify_params)

        #must add raster to iface so that is becomes active layer, then symbolize it in next step
        output_raster = QgsRasterLayer(output_raster_destination)
        iface.addRasterLayer(output_raster_destination)

        #grabs active layer and data from that layer
        layer = iface.activeLayer()
        provider = layer.dataProvider()
        extent = layer.extent()
        
        raster_stats = provider.bandStatistics(1, QgsRasterBandStats.All) 




        #   // STEP 4. Symbolize output layer //
        log("Symbolizing Output Layer")

        #creates raster shader and creates discrete color ramp
        raster_shader = QgsColorRampShader()
        raster_shader.setColorRampType(QgsColorRampShader.Discrete)           

        #creates layer symbology from raster_stats data
        raster_stats = provider.bandStatistics(1, QgsRasterBandStats.All)
        symbology = Symbology(raster_stats, input_esv_field)
        colors_list = symbology.symbolize_input_raster()

        raster_shader.setColorRampItemList(colors_list)         #applies symbology to raster_shader
        shader = QgsRasterShader()
        shader.setRasterShaderFunction(raster_shader)       

        renderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(), 1, shader)    #renders selected raster layer
        layer.setRenderer(renderer)
        layer.triggerRepaint()
        

    
        log(self.tr(f"Adding final raster to map."))
        #need to add result from gdal:rastercalculator to map (doesn't happen automatically)
        
        
        log(self.tr("Done!\n"))



        # 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 result


    def flags(self):
        """
        From documentation: Algorithm is not thread safe and cannot be run in a
        background thread, e.g. algorithms which manipulate the current project,
        layer selections, or with external dependencies which are not thread safe.
        """
        return super().flags() | QgsProcessingAlgorithm.FlagNoThreading

    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 'Step 2: Map the value of individual ecosystem services'

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

    def shortHelpString(self):
        """
        Returns a localised short helper string for the algorithm. This string
        should provide a basic description about what the algorithm does and the
        parameters and outputs associated with it..
        """
        return self.tr("This algorithm takes the clipped Land Use/Land Cover (LULC) raster and Input ESV table from step 1 and creates a new raster for which the value is the corresponding per-pixel value (minimum, mean, or maximum) of the user-chosen ecosystem service. The new raster is then colored according to the ecosystem service chosen. It's values are divided into even quintiles to emphasize breaks in the data. The user can repeat this step for additional levels (min, mean, max) and ecosystem services.\n ~~~~~~~~~~~~~~~~ \n Inputs: \n Select land use/land cover data source: Choose either NLCD (National Land Cover Dataset) or NALCMS (North American Land Change Monitoring System). \n Input clipped raster layer: This should be your output from step 1, which is a LULC raster layer clipped to the extent of your study area. \n Choose ecosystem service of interest: Specify the ecosystem service you are interested in. \n Choose ecosystem Service Value Level: Choose to map minimum, mean, or maximum values from the ESV table for your corresponding ESV choice. \n Output esv Raster: Specify an output location for your ESV raster. THIS CANNOT BE BLANK \n See “Help” for more information on value origins and ecosystem service descriptions.")

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

    def helpUrl(self):
        return "http://www.keylogeconomics.com/ecovaluator.html"

    def createInstance(self):
        return MapTheValueOfIndividualEcosystemServices()