CostRasterCreator_algorithm.py

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

"""
/***************************************************************************
 ForestRoads
                                 A QGIS plugin
 Create a network of forest roads based on zones to access, roads to connect
 them to, and a cost matrix.
 The code of the plugin is based on the "LeastCostPath" plugin available on
 https://github.com/Gooong/LeastCostPath. We thank their team for the template.
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 10-07-2019
        copyright            : (C) 2019 by Clement Hardy
        email                : clem.hardy@outlook.fr
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
 This script describes the algorithm used to help the user to create the cost raster that
 will be used for the forest road network algorithm.
"""

__author__ = 'clem.hardy@outlook.fr'
__date__ = 'Currently in work'
__copyright__ = '(C) 2019 by Clement Hardy'

# We load every function necessary from the QIS packages.
import os
import random
import queue
import math
from PyQt5.QtCore import QCoreApplication, QVariant
from PyQt5.QtGui import QIcon
from qgis.core import (
    Qgis,
    QgsRasterBlock,
    QgsRasterFileWriter,
    QgsFeature,
    QgsGeometry,
    QgsPoint,
    QgsPointXY,
    QgsField,
    QgsFields,
    QgsWkbTypes,
    QgsProcessing,
    QgsFeatureSink,
    QgsProcessingException,
    QgsProcessingAlgorithm,
    QgsProcessingParameterRasterDestination,
    QgsProcessingParameterField,
    QgsProcessingParameterFeatureSink,
    QgsProcessingParameterRasterLayer,
    QgsProcessingParameterBand,
    QgsProcessingParameterBoolean,
    QgsProcessingParameterNumber,
    QgsProcessingParameterEnum,
    QgsProcessingParameterMatrix,
    QgsRasterFileWriter
)
# We import mathematical functions needed for the algorithm.
from math import floor, sqrt


# The algorithm class heritates from the algorithm class of QGIS.
# There, it can register different parameter during initialization
# that can be put into variables using "
class CostRasterAlgorithm(QgsProcessingAlgorithm):
    """
    Class that described the algorithm, that will be launched
    via the provider, itself launched via initialization of
    the plugin.
    """

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

    INITIAL_ROAD_NETWORK = 'INITIAL_ROAD_NETWORK'

    BASIC_DISTANCE_COST = 'BASIC_DISTANCE_COST'

    COARSE_ELEVATION_RASTER = 'COARSE_ELEVATION_RASTER'

    COARSE_ELEVATION_COSTS = 'COARSE_ELEVATION_COSTS'

    FINE_ELEVATION_RASTER = 'FINE_ELEVATION_RASTER'

    FINE_ELEVATION_COSTS = 'FINE_ELEVATION_COSTS'

    COARSE_WATER_RASTER = 'COARSE_WATER_RASTER'

    COARSE_WATER_COST = 'COARSE_WATER_COST'

    FINE_WATER_RASTER = 'FINE_WATER_RASTER'

    FINE_WATER_COST = 'FINE_WATER_COST'

    SOIL_RASTER = 'SOIL_RASTER'

    ADDITIONAL_COST_RASTER = 'ADDITIONAL_COST_RASTER'

    OUTPUT = 'OUTPUT'

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

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INITIAL_ROAD_NETWORK,
                self.tr('Initial road network raster'),
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.BASIC_DISTANCE_COST,
                self.tr(
                    'Basic distance cost'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=10000,
                optional=False,
                minValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.COARSE_ELEVATION_RASTER,
                self.tr('Coarse elevation raster'),
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterMatrix(self.COARSE_ELEVATION_COSTS,
                                         "Additional costs associated with mean values of slope in a pixel",
                                         numberRows=2,
                                         hasFixedNumberRows=False,
                                         headers=["Lower slope threshold", "Upper slope threshold", "Additional Cost"],
                                         defaultValue = ['0', '10', '0', '10', '100', '5000'],
                                         optional = True)
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.FINE_ELEVATION_RASTER,
                self.tr('Fine elevation raster'),
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterMatrix(self.FINE_ELEVATION_COSTS,
                                         "Multiplicative costs associated with values of fine topography in a pixel",
                                         numberRows=2,
                                         hasFixedNumberRows=False,
                                         headers=["Lower threshold", "Upper threshold", "Multiplicative Cost"],
                                         defaultValue=['0', '300', '1', '300', '1000', '2'],
                                         optional = True)
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.COARSE_WATER_RASTER,
                self.tr('Coarse Water raster'),
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.COARSE_WATER_COST,
                self.tr(
                    'Coarse water cost'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=500000,
                optional=True,
                minValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.FINE_WATER_RASTER,
                self.tr('Fine Water raster'),
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.FINE_WATER_COST,
                self.tr(
                    'Fine water cost'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=20000,
                optional=True,
                minValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.SOIL_RASTER,
                self.tr('Soil raster'),
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.ADDITIONAL_COST_RASTER,
                self.tr('Additional cost raster to avoid certain areas (e.g. sensitive flora or fauna zones, etc.)'),
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterRasterDestination(
                self.OUTPUT,
                self.tr('Output file of the algorithm')
            )
        )

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

        if self.parameterAsRasterLayer(
            parameters,
            self.INITIAL_ROAD_NETWORK,
            context
        ):
            initial_road_network_raster = self.parameterAsRasterLayer(
                parameters,
                self.INITIAL_ROAD_NETWORK,
                context
            )
        else:
            initial_road_network_raster = None

        basic_distance_cost = self.parameterAsDouble(
            parameters,
            self.BASIC_DISTANCE_COST,
            context
        )

        if self.parameterAsRasterLayer(
            parameters,
            self.COARSE_ELEVATION_RASTER,
            context
        ):
            coarse_elevation_raster = self.parameterAsRasterLayer(
                parameters,
                self.COARSE_ELEVATION_RASTER,
                context
            )
        else:
            coarse_elevation_raster = None

        if self.parameterAsMatrix(
            parameters,
            self.COARSE_ELEVATION_COSTS,
            context
        ):
            coarse_elevation_costs = self.parameterAsMatrix(
                parameters,
                self.COARSE_ELEVATION_COSTS,
                context
            )
            coarse_elevation_costs = CostRasterCreatorHelper.CollapsedTableToMatrix(coarse_elevation_costs, 3)
        else:
            coarse_elevation_costs = None

        if self.parameterAsRasterLayer(
            parameters,
            self.FINE_ELEVATION_RASTER,
            context
        ):
            fine_elevation_raster = self.parameterAsRasterLayer(
                parameters,
                self.FINE_ELEVATION_RASTER,
                context
            )
        else:
            fine_elevation_raster = None

        if self.parameterAsMatrix(
            parameters,
            self.FINE_ELEVATION_COSTS,
            context
        ):
            fine_elevation_costs = self.parameterAsMatrix(
                parameters,
                self.FINE_ELEVATION_COSTS,
                context
            )
            fine_elevation_costs = CostRasterCreatorHelper.CollapsedTableToMatrix(fine_elevation_costs, 3)
        else:
            fine_elevation_costs = None

        if self.parameterAsRasterLayer(
            parameters,
            self.COARSE_WATER_RASTER,
            context
        ):
            coarse_water_raster = self.parameterAsRasterLayer(
                parameters,
                self.COARSE_WATER_RASTER,
                context
            )
        else:
            coarse_water_raster = None

        if self.parameterAsDouble(
            parameters,
            self.COARSE_WATER_COST,
            context
        ):
            coarse_water_cost = self.parameterAsDouble(
                parameters,
                self.COARSE_WATER_COST,
                context
            )
        else:
            coarse_water_cost = None

        if self.parameterAsRasterLayer(
            parameters,
            self.FINE_WATER_RASTER,
            context
        ):
            fine_water_raster = self.parameterAsRasterLayer(
                parameters,
                self.FINE_WATER_RASTER,
                context
            )
        else:
            fine_water_raster = None

        if self.parameterAsDouble(
            parameters,
            self.FINE_WATER_COST,
            context
        ):
            fine_water_cost = self.parameterAsDouble(
                parameters,
                self.FINE_WATER_COST,
                context
            )
        else:
            fine_water_cost = None

        if self.parameterAsRasterLayer(
            parameters,
            self.SOIL_RASTER,
            context
        ):
            soil_raster = self.parameterAsRasterLayer(
                parameters,
                self.SOIL_RASTER,
                context
            )
        else:
            soil_raster = None

        if self.parameterAsRasterLayer(
            parameters,
            self.ADDITIONAL_COST_RASTER,
            context
        ):
            special_raster = self.parameterAsRasterLayer(
                parameters,
                self.ADDITIONAL_COST_RASTER,
                context
            )
        else:
            special_raster = None

        feedback.pushInfo(self.tr("Checking inputs..."))
        # If source was not found, throw an exception to indicate that the algorithm
        # encountered a fatal error. The exception text can be any string, but in this
        # case we use the pre-built invalidSourceError method to return a standard
        # helper text for when a source cannot be evaluated
        if basic_distance_cost is None or basic_distance_cost == 0:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.BASIC_DISTANCE_COST))

        # Now, we check to see if the CRS, extent and resolution of every raster is equal, so that their align properly.
        listOfRastersToCheck = list()
        if initial_road_network_raster is not None:
            listOfRastersToCheck.append(initial_road_network_raster)
        if coarse_elevation_raster is not None:
            listOfRastersToCheck.append(coarse_elevation_raster)
        if fine_elevation_raster is not None:
            listOfRastersToCheck.append(fine_elevation_raster)
        if coarse_water_raster is not None:
            listOfRastersToCheck.append(coarse_water_raster)
        if fine_water_raster is not None:
            listOfRastersToCheck.append(fine_water_raster)
        if soil_raster is not None:
            listOfRastersToCheck.append(soil_raster)
        if special_raster is not None:
            listOfRastersToCheck.append(special_raster)

        if len(listOfRastersToCheck) == 0:
            raise QgsProcessingException(self.tr("At least one input raster is needed ! Please, input one raster."))

        # We check that every raster has the same CRS, extent and resolution.
        CostRasterCreatorHelper.CheckRastersCompatibility(listOfRastersToCheck)

        # We also check that the matrix of parameters entered for the coarse elevation costs and fine elevation costs
        # are correct : meaning that there are no holes between the thresholds, and that every lower threshold is lower
        # than the upper threshold.
        if coarse_elevation_costs is not None:
            CostRasterCreatorHelper.CheckThresholds(coarse_elevation_costs, "Coarse elevation costs")
        if fine_elevation_costs is not None:
            CostRasterCreatorHelper.CheckThresholds(fine_elevation_costs, "Fine elevation costs")

        # Then, we create the raster blocks
        if initial_road_network_raster is not None:
            initial_road_network_raster_block = CostRasterCreatorHelper.get_all_block(initial_road_network_raster)
        else:
            initial_road_network_raster_block = None

        if coarse_elevation_raster is not None:
            coarse_elevation_raster_block = CostRasterCreatorHelper.get_all_block(coarse_elevation_raster)
        else:
            coarse_elevation_raster_block = None

        if fine_elevation_raster is not None:
            fine_elevation_raster_block = CostRasterCreatorHelper.get_all_block(fine_elevation_raster)
        else:
            fine_elevation_raster_block = None

        if coarse_water_raster is not None:
            coarse_water_raster_block = CostRasterCreatorHelper.get_all_block(coarse_water_raster)
        else:
            coarse_water_raster_block = None

        if fine_water_raster is not None:
            fine_water_raster_block = CostRasterCreatorHelper.get_all_block(fine_water_raster)
        else:
            fine_water_raster_block = None

        if soil_raster is not None:
            soil_raster_block = CostRasterCreatorHelper.get_all_block(soil_raster)
        else:
            soil_raster_block = None

        if special_raster is not None:
            special_raster_block = CostRasterCreatorHelper.get_all_block(special_raster)
        else:
            special_raster_block = None

        feedback.pushInfo(self.tr("Preparing output..."))
        # We set the output to be ready : It is a QgsDataProvider
        outputFile = self.parameterAsOutputLayer(parameters, self.OUTPUT, context)
        outputFormat = QgsRasterFileWriter.driverForExtension(os.path.splitext(outputFile)[1])

        crs = listOfRastersToCheck[0].crs()
        extent = listOfRastersToCheck[0].extent()
        rows = max([math.ceil(extent.height() / listOfRastersToCheck[0].rasterUnitsPerPixelY()), 1.0])
        cols = max([math.ceil(extent.width() / listOfRastersToCheck[0].rasterUnitsPerPixelX()), 1.0])
        # We will need this value for the fine water computation later
        pixelSide = (listOfRastersToCheck[0].rasterUnitsPerPixelY() +
                             listOfRastersToCheck[0].rasterUnitsPerPixelX()) / 2

        writer = QgsRasterFileWriter(outputFile)
        writer.setOutputProviderKey('gdal')
        writer.setOutputFormat(outputFormat)
        provider = writer.createOneBandRaster(Qgis.Float32, cols, rows, extent, crs)

        if provider is None:
            raise QgsProcessingException(self.tr("Could not create raster output: {}").format(outputFile))
        if not provider.isValid():
            raise QgsProcessingException(self.tr("Could not create raster output {}").format(outputFile))

        provider.setNoDataValue(1, -9999)

        # We create the data block for the output raster, and we fill it with the values we need
        dataBlock = QgsRasterBlock(Qgis.Float32, cols, rows)
        # i = float(1.0)

        feedback.pushInfo(self.tr("Calculating cost raster..."))
        progress = 0
        feedback.setProgress(0)
        errorMessages = list()
        for y in range(dataBlock.height()):
            for x in range(dataBlock.width()):
                if feedback.isCanceled():
                    raise QgsProcessingException(self.tr("ERROR: Operation was cancelled."))

                    # If there is a road already on this pixel, then the cost is 0.
                if initial_road_network_raster_block is not None and \
                     (initial_road_network_raster_block.value(y,
                                                              x) != 0 and not initial_road_network_raster_block.isNoData(
                         y, x)):
                    finalValue = 0

                # If there is a water body on this pixel, we stop everything :
                # The cost will be the construction of a bridge
                elif coarse_water_raster_block is not None and coarse_water_raster_block.value(y, x) != 0:
                    if coarse_water_cost is not None:
                        finalValue = coarse_water_cost
                    else:
                        raise QgsProcessingException(self.tr("A coarse water raster has been given, but no coarse " +
                                                             "water cost. Please input a coarse water cost."))
                    # feedback.pushInfo("Seems like there was a road on this pixel. Final value is " + str(finalValue))

                # Else, if we are not on a road or on a body of water...
                else:
                    # We start with the base cost of crossing the pixel
                    finalValue = basic_distance_cost
                    # feedback.pushInfo("No road on this pixel, we put basic distance cost. Final value is " + str(finalValue))
                    # Then, if we have it, we add the soil cost
                    if soil_raster_block is not None:
                        finalValue += soil_raster_block.value(y, x)
                        # feedback.pushInfo("After soils, final value is " + str(finalValue))
                    if special_raster_block is not None:
                        finalValue += special_raster_block.value(y, x)
                    # Then the coarse elevation value
                    if coarse_elevation_raster_block is not None:
                        additionalValue, errorMessage = CostRasterCreatorHelper.CalculateCoarseElevationCost(y,
                                                                                           x,
                                                                                           coarse_elevation_raster_block,
                                                                                           coarse_elevation_costs,
                                                                                           pixelSide)
                        finalValue += additionalValue
                        if errorMessage is not None:
                            errorMessages.append(errorMessage)
                        # feedback.pushInfo("After coarse elevation, final value is " + str(finalValue))
                    # Then the fine water value
                    if fine_water_raster_block is not None:
                        finalValue += CostRasterCreatorHelper.CalculateFineWaterCost(y,
                                                                                     x,
                                                                                     fine_water_raster_block,
                                                                                     fine_water_cost,
                                                                                     pixelSide)
                        # feedback.pushInfo("After fine water, final value is " + str(finalValue))
                    # Then, we multiply everything with the fine elevation cost.
                    if fine_elevation_raster_block is not None:
                        finalValue, errorMessage = CostRasterCreatorHelper.CalculateFineElevationCost(y,
                                                                                    x,
                                                                                    fine_elevation_raster_block,
                                                                                    fine_elevation_costs,
                                                                                    finalValue)
                        if errorMessage is not None:
                            errorMessages.append(errorMessage)
                        # feedback.pushInfo("After fine elevation, final value is " + str(finalValue))
                dataBlock.setValue(y, x, float(finalValue))

                progress += 1
                feedback.setProgress(100 * (progress / (dataBlock.height()*dataBlock.width())))

        # We write the values in the provider of our files
        provider.writeBlock(dataBlock, 1)

        # We stop the edition of the output raster
        provider.setEditable(False)

        # We display the warning messages about the thresholds
        if len(errorMessages) > 0:
            above = 0
            below = 0
            for errorMessage in errorMessages:
                if errorMessage == "Above":
                    above += 1
                elif errorMessage == "Below":
                    below += 1
            feedback.pushInfo(self.tr("WARNING : There were " + str(below) + " situations where the value of a pixel was under"
                            " the lowest threshold given as a parameter; and " + str(above) + " when it was above the"
                            " highest. In those cases, the lowest or highest value of the parameter range was used."
                            " To avoid that, please make sure to use thresholds that cover all of the range of values in"
                            " your rasters."))


        # We make the output
        return {self.OUTPUT: outputFile}

    # Here are different functions used by QGIS to name and define the algorithm
    # to the user.
    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 'Cost Raster Creator'

    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.
        """
        # Not enough algorithms in this plugin for the need to make groups of algorithms
        return ''

    # Function used for translation. Called every time something needs to be
    # Translated
    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return CostRasterAlgorithm()

    def helpUrl(self):
        return 'https://github.com/Klemet/ForestRoadNetworkPluginForQGIS'

    def shortHelpString(self):
        return self.tr("""
        This algorithm allows you to create a cost raster for the "Forest Road Network Creation" algorithm using several rasters and cost values. Each pixel of the raster will represent the cost of crossing the said raster. At least one raster/associated cost value and the base distance cost are needed; all of the rest is optional. However, the more rasters and costs, the more realistic the resulting cost raster will be. The output will contain all of the input cost and rasters in its metadata to allow the user to easily retrieve information about its creation.

        WARNING : The algorithm assumes that the raster's data is contained in the first raster band of the rasters.
        
        **Parameters:**

          Please ensure all the input layers have the same CRS, extent and resolution; the output raster will share those characteristrics.

          - Initial road network raster : This raster should contain information about the presence or absence of roads in its pixels. 0 = no roads, any other value = presence of a road. In the resulting cost raster, passing by a pixel with an already existing road will be free. We recommand the use of the "rasterize" QGIS tool with a vector layer of existing roads to create this raster.

          - Basic distance cost : The irreducible basic cost of constructing a forest road on the distance of the side of a pixel. We recommand that it is expressed as a monetary unit (dollards, etc.)

          - Coarse elevation raster : A raster containing the mean elevation of each pixel. Basically a DEM (Digital Elevation Model)

          - Coarse elevation cost : A value that will be multiplied with the percentage of slope between the neighbouring pixels; the resulting number will be added to the cost of crossing this pixel (addition). Represents the augmentation of the cost to construct a road on a slope.

          - Fine elevation raster : A raster containing the number of contour/elevation lines, or their added length inside each pixel. We recommand the use of the combine use of the "Create Grid", "Sum Line Length" and "Rasterize" tools in QGIS to create it from a file containing contour/elevation lines.

          - Fine elevation threshold : The threshold of fine elevation value above which the contour lines are so numerous in the pixel that it is considered that detours will have to be made, multiplying the distance of road constructed in this pixel by 2. A proportional relation is applied for values fine elevation values higher than this threshold.

          - Coarse water raster : A raster containing either 1 if the pixel is part of a big body of water (river or lake); or 0 if that is not the case.
          
          - Coarse water cost : The cost of building a bridge on the distance of one pixel on a big body of water.
          
          - Fine water raster : A raster containing the added length of the stream crossing the raster. To create it, we recommand the same approach as for the Fine Elevation Raster. The number of culverts that will have to be created if crossing the pixel will be deduced by the algorithm.
          
          - Fine water cost : The cost of building a culvert. It will be added to the cost of building a road on the given pixel.
          
          - Soil raster : A raster directly containing an additional cost of building a road on this pixel due to the soil of the pixel. It can represent the added cost of having to bring gravel from afar.
        
          - Additional cost raster : A raster containing additional cost to avoid certain areas in particular. Examples can be particular fauna habitats, the edge of lakes, etc.        
        """)

    def shortDescription(self):
        return self.tr('Create a cost raster representing the cost of constructing a forest road on a given pixel.')

    # Path to the icon of the algorithm
    def svgIconPath(self):
        return '.icon.png'

    def tags(self):
        return ['wood', 'flux', 'type', 'primary', 'secondary', 'roads', 'analysis', 'road', 'network',
                'forest', 'tertiary', 'temporary']

# Methods to help the algorithm; all static, do not need to initialize an object of this class.
class CostRasterCreatorHelper:

    @staticmethod
    def CollapsedTableToMatrix(matrixOfThresholds, numberOfColumns):
        """This function exist to counter the fact that QGIS collapses the matrix of parameters into a list. It puts
        it back into a proper matrix, meaning a list of lists. It also transforms the strings in the matrix into numbers."""

        matrixOfThresholdsAsNumber = list(map(float, matrixOfThresholds))

        numberOfRows = len(matrixOfThresholdsAsNumber) // numberOfColumns
        matrixToReturn = list()

        for indexing in range(numberOfRows):
            matrixToReturn.append(matrixOfThresholdsAsNumber[(indexing * numberOfColumns):(indexing * numberOfColumns + numberOfColumns)])

        return matrixToReturn

    # Function that get the data block from a entire raster for a given band
    @staticmethod
    def get_all_block(raster_layer):
        provider = raster_layer.dataProvider()
        extent = provider.extent()

        xres = raster_layer.rasterUnitsPerPixelX()
        yres = raster_layer.rasterUnitsPerPixelY()
        width = floor((extent.xMaximum() - extent.xMinimum()) / xres)
        height = floor((extent.yMaximum() - extent.yMinimum()) / yres)
        return provider.block(1, extent, width, height)

    @staticmethod
    def CheckRastersCompatibility(listOfRasterLayers):
        """This function reads the characteristics of several rasters in a list,
        and raise an exception if the CRS, extent or resolution are not the same for each of them."""

        for raster_layer in listOfRasterLayers:
            for another_raster_layer in listOfRasterLayers:
                if raster_layer != another_raster_layer:
                    if raster_layer.crs() != another_raster_layer.crs():
                        raise QgsProcessingException("ERROR: The input rasters have different CRSs.")
                    if raster_layer.extent() != another_raster_layer.extent():
                        raise QgsProcessingException("ERROR: The input rasters have different extents.")
                    if raster_layer.height() != another_raster_layer.height():
                        raise QgsProcessingException("ERROR: The input rasters have different heights.")
                    if raster_layer.width() != another_raster_layer.width():
                        raise QgsProcessingException("ERROR: The input rasters have different widths.")
                    if raster_layer.rasterUnitsPerPixelX() != another_raster_layer.rasterUnitsPerPixelX():
                        raise QgsProcessingException("ERROR: The input rasters have different horizontal resolution.")
                    if raster_layer.rasterUnitsPerPixelY() != another_raster_layer.rasterUnitsPerPixelY():
                        raise QgsProcessingException("ERROR: The input rasters have different vertical resolution.")

    @staticmethod
    def CheckThresholds(matrixOfThresholds, thresholdsName):
        """This function check if there are no holes in threshold parameters, and if every lower threshold is inferior
        to the associated upper threshold. If not, we raise an exception."""

        for index in range(len(matrixOfThresholds)):
            lowerThreshold = matrixOfThresholds[index][0]
            upperThreshold = matrixOfThresholds[index][1]
            if lowerThreshold >= upperThreshold:
                raise QgsProcessingException("ERROR: For the parameter " + str(thresholdsName) + ", the lower threshold of"
                                            " row " + str(index) + " is bigger or equal to the upper threshold of the same row."
                                            " Please, correct it before trying again.")
            if (index + 1) < len(matrixOfThresholds):
                lowerThresholdOfHigherRow = matrixOfThresholds[index + 1][0]
                if (lowerThresholdOfHigherRow - upperThreshold) != 0:
                    raise QgsProcessingException(
                        "ERROR: For the parameter " + str(thresholdsName) + ", there is a hole between the upper threshold of"
                        " row " + str(index) + " and the lower threshold of row " + str(index+1) + ". The thresholds must"
                        " be continuous and ordered from smaller to higher. Please, correct it before trying again.")

    @staticmethod
    def GetThresholdValue(matrixOfThresholds, valueInsideThresholds):
        """This function gets the associated value (additional value or multiplicative value) inside a matrix of thresholds
        (third column) given a value that will be inside one of the thresholds. Is also register error messages to warn
        the user about values that are not into its thresholds."""

        thresholdValue = None
        errorMessage = None

        # First, we treat the case of our value being under the lowest threshold, or higher than the highest threshold.
        if valueInsideThresholds <  matrixOfThresholds[0][0]:
            thresholdValue = matrixOfThresholds[0][2]
            errorMessage = "Below"
        elif valueInsideThresholds > matrixOfThresholds[len(matrixOfThresholds) - 1][1]:
            thresholdValue = matrixOfThresholds[len(matrixOfThresholds) - 1][2]
            errorMessage = "Above"
        else:
            for index in range(len(matrixOfThresholds)):
                if valueInsideThresholds >= matrixOfThresholds[index][0] and valueInsideThresholds < matrixOfThresholds[index][1]:
                    thresholdValue = matrixOfThresholds[index][2]
                    break

        if thresholdValue is None:
            raise QgsProcessingException("ERROR: Couldn't find value " + str(valueInsideThresholds) + " in the following"
                                         " matrix of thresholds : " + str(matrixOfThresholds))
        else:
            return thresholdValue, errorMessage

    @staticmethod
    def CalculateCoarseElevationCost(row, column, coarse_elevation_raster_block, coarse_elevation_costs, pixelSide):
        """"A function to calculate the coarse elevation cost of a pixel. Current method is to make the mean of the
        difference in elevation between neighbouring pixels, to transform it into slope percentage using the pixel size,
        and then multiply it by the coarse elevation cost."""

        meanSlope = 0
        numberOfNeighbors = 0
        cellElevation = coarse_elevation_raster_block.value(row, column)

        # North-west neighbor
        if row - 1 >= 0 and column - 1 >= 0:
            elevationDifference = abs(coarse_elevation_raster_block.value(row - 1, column - 1) - cellElevation)
            horizontalDistance = math.sqrt(2) * pixelSide
            meanSlope += (elevationDifference/horizontalDistance) * 100
            numberOfNeighbors += 1
        # North neighbor
        if row - 1 >= 0:
            elevationDifference = abs(coarse_elevation_raster_block.value(row - 1, column) - cellElevation)
            horizontalDistance = 1 * pixelSide
            meanSlope += (elevationDifference / horizontalDistance) * 100
            numberOfNeighbors += 1
        # North-east neighbour
        if row - 1 >= 0 and column + 1 < coarse_elevation_raster_block.width():
            elevationDifference = abs(coarse_elevation_raster_block.value(row - 1, column + 1) - cellElevation)
            horizontalDistance = math.sqrt(2) * pixelSide
            meanSlope += (elevationDifference / horizontalDistance) * 100
            numberOfNeighbors += 1
        # Western neighbor
        if column - 1 >= 0:
            elevationDifference = abs(coarse_elevation_raster_block.value(row, column - 1) - cellElevation)
            horizontalDistance = 1 * pixelSide
            meanSlope += (elevationDifference / horizontalDistance) * 100
            numberOfNeighbors += 1
        # Eastern neighbor
        if column + 1 < coarse_elevation_raster_block.width():
            elevationDifference = abs(coarse_elevation_raster_block.value(row, column + 1) - cellElevation)
            horizontalDistance = 1 * pixelSide
            meanSlope += (elevationDifference / horizontalDistance) * 100
            numberOfNeighbors += 1
        # Southern neighbor
        if row + 1 < coarse_elevation_raster_block.height():
            elevationDifference = abs(coarse_elevation_raster_block.value(row + 1, column) - cellElevation)
            horizontalDistance = 1 * pixelSide
            meanSlope += (elevationDifference / horizontalDistance) * 100
            numberOfNeighbors += 1
        # South-west neighbor
        if row + 1 < coarse_elevation_raster_block.height() and column - 1 >= 0:
            elevationDifference = abs(coarse_elevation_raster_block.value(row + 1, column - 1) - cellElevation)
            horizontalDistance = math.sqrt(2) * pixelSide
            meanSlope += (elevationDifference / horizontalDistance) * 100
            numberOfNeighbors += 1
        # South-east neighbor
        if row + 1 < coarse_elevation_raster_block.height() and column + 1 < coarse_elevation_raster_block.width():
            elevationDifference = abs(coarse_elevation_raster_block.value(row + 1, column + 1) - cellElevation)
            horizontalDistance = math.sqrt(2) * pixelSide
            meanSlope += (elevationDifference / horizontalDistance) * 100
            numberOfNeighbors += 1

        meanSlope = meanSlope / float(numberOfNeighbors)

        additionalCost, errorMessage = CostRasterCreatorHelper.GetThresholdValue(coarse_elevation_costs, meanSlope)
        return meanSlope * additionalCost, errorMessage

    @staticmethod
    def CalculateFineWaterCost(row, column, fine_water_raster_block, fine_water_cost, pixelSize):
        """"A function to calculate the fine water cost of a pixel. For that, we take the size of a pixel, and
        consider how many streams will have to be crossed if their are some in the cell."""

        distanceOfStreamInPixel = fine_water_raster_block.value(row, column)
        numberOfStreamsToCross = distanceOfStreamInPixel // pixelSize

        return numberOfStreamsToCross * fine_water_cost

    @staticmethod
    def CalculateFineElevationCost(row, column, fine_elevation_raster_block, fine_elevation_costs, valueToUpdate):
        """A function to calculate the fine elevation cost of a pixel. We will take the cost for this pixel, and
        multiply it by how much detour must be done because of fine changes in the topography."""

        # The threshold given by the user is for multiplying the cost of construction by 2. We would it for 1, so as
        # to make a floor division in the case of higher values.
        # singleUnitOfRoadThreshold = fine_elevation_threshold / float(2)

        #multiplicationDueToFineElevation = fine_elevation_raster_block.value(row, column) // singleUnitOfRoadThreshold
        fineElevationValue = fine_elevation_raster_block.value(row, column)
        multiplicationDueToFineElevation, errorMessage = CostRasterCreatorHelper.GetThresholdValue(fine_elevation_costs, fineElevationValue)

        if multiplicationDueToFineElevation > 1:
            return valueToUpdate * multiplicationDueToFineElevation, errorMessage
        else:
            return valueToUpdate, errorMessage