forestRoadNetwork_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 make the forest road network.
"""

__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 random
from .kdtree import KDTree
import numpy as np
from PyQt5.QtCore import QCoreApplication, QVariant
from PyQt5.QtGui import QIcon
from qgis.core import (
    QgsFeature,
    QgsGeometry,
    QgsPoint,
    QgsPointXY,
    QgsField,
    QgsFields,
    QgsWkbTypes,
    QgsProcessing,
    QgsFeatureSink,
    QgsProcessingException,
    QgsProcessingAlgorithm,
    QgsProcessingParameterFeatureSource,
    QgsProcessingParameterFeatureSink,
    QgsProcessingParameterRasterLayer,
    QgsProcessingParameterField,
    QgsProcessingParameterBand,
    QgsProcessingParameterBoolean,
    QgsProcessingParameterNumber,
    QgsProcessingParameterEnum
)
# We import the algorithm used for processing a road.
from .dijkstra_algorithm import dijkstra
# 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 ForestRoadNetworkAlgorithm(QgsProcessingAlgorithm):
    """
    Class that described the algorithm, that will be launched
    via the provider, itself launched via initialization of
    the plugin.

    The algorithm takes 4 entries :

    - A cost raster
    - The raster band to use for the cost
    - The layer with the polygons of zones to access
    - The layer with the roads (lines) that they can be connected to
    by the generated roads
    """

    # 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_COST_RASTER = 'INPUT_COST_RASTER'

    INPUT_RASTER_BAND = 'INPUT_RASTER_BAND'

    INPUT_POLYGONS_TO_ACCESS = 'INPUT_POLYGONS_TO_ACCESS'

    INPUT_ROADS_TO_CONNECT_TO = 'INPUT_ROADS_TO_CONNECT_TO'

    SKIDDING_DISTANCE = 'SKIDDING_DISTANCE'

    METHOD_OF_GENERATION = 'METHOD_OF_GENERATION'

    HEURISTIC_IN_POLYGONS = 'HEURISTIC_IN_POLYGONS'

    ANGLES_CONSIDERED = 'ANGLES_CONSIDERED'

    PUNISHER_45DEGREES = 'PUNISHER_45DEGREES'

    PUNISHER_90DEGREES = 'PUNISHER_90DEGREES'

    PUNISHER_135DEGREES = 'PUNISHER_135DEGREES'

    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.INPUT_COST_RASTER,
                self.tr('Cost raster layer'),
            )
        )

        self.addParameter(
            QgsProcessingParameterBand(
                self.INPUT_RASTER_BAND,
                self.tr('Cost raster band'),
                0,
                self.INPUT_COST_RASTER,
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT_POLYGONS_TO_ACCESS,
                self.tr('Polygons to access via the generated roads'),
                [QgsProcessing.TypeVectorPolygon]
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT_ROADS_TO_CONNECT_TO,
                self.tr('Roads to connect the polygons to access to'),
                [QgsProcessing.TypeVectorLine]
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.SKIDDING_DISTANCE,
                self.tr('Skidding distance (in CRS units)'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=100,
                optional=False,
                minValue=0
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.METHOD_OF_GENERATION,
                self.tr('Method of generation of the road network'),
                ['Random', 'Closest first', 'Farthest first'],
                defaultValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.HEURISTIC_IN_POLYGONS,
                self.tr('Attribute field of the harvest polygons that indicates the order in which they must be reached'),
                parentLayerParameterName=self.INPUT_POLYGONS_TO_ACCESS,
                type=QgsProcessingParameterField.Numeric,
                optional=True
            )
        )

        self.addParameter(
            QgsProcessingParameterEnum(
                self.ANGLES_CONSIDERED,
                self.tr('Should the angles of the roads be considered in the cost ?'),
                ['Yes', 'No'],
                defaultValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.PUNISHER_45DEGREES,
                self.tr('Punishing multiplier for a 45 degrees turning angle'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=1.25,
                optional=True,
                minValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.PUNISHER_90DEGREES,
                self.tr('Punishing multiplier for a 90 degrees turning angle'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=2,
                optional=True,
                minValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.PUNISHER_135DEGREES,
                self.tr('Punishing multiplier for a 135 degrees turning angle'),
                type=QgsProcessingParameterNumber.Double,
                defaultValue=5,
                optional=True,
                minValue=1
            )
        )

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

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

        cost_raster_band = self.parameterAsInt(
            parameters,
            self.INPUT_RASTER_BAND,
            context
        )

        polygons_to_connect = self.parameterAsVectorLayer(
            parameters,
            self.INPUT_POLYGONS_TO_ACCESS,
            context
        )

        current_roads = self.parameterAsVectorLayer(
            parameters,
            self.INPUT_ROADS_TO_CONNECT_TO,
            context
        )

        skidding_distance = self.parameterAsInt(
            parameters,
            self.SKIDDING_DISTANCE,
            context
        )

        method_of_generation = self.parameterAsString(
            parameters,
            self.METHOD_OF_GENERATION,
            context
        )

        if self.parameterAsString(
                    parameters,
                    self.HEURISTIC_IN_POLYGONS,
                    context
                ) is not None:
            heuristic_in_polygons_index = polygons_to_connect.fields().lookupField(
                self.parameterAsString(
                    parameters,
                    self.HEURISTIC_IN_POLYGONS,
                    context
                ))
        else:
            heuristic_in_polygons_index = None

        if self.parameterAsString(
            parameters,
            self.ANGLES_CONSIDERED,
            context
        ) == '0' :
            angles_considered = True
        else:
            angles_considered = False

        punisherAngleDictionnary = dict()

        punisherAngleDictionnary[45] = self.parameterAsDouble(
            parameters,
            self.PUNISHER_45DEGREES,
            context
        )

        punisherAngleDictionnary[90] = self.parameterAsDouble(
            parameters,
            self.PUNISHER_90DEGREES,
            context
        )

        punisherAngleDictionnary[135] = self.parameterAsDouble(
            parameters,
            self.PUNISHER_135DEGREES,
            context
        )

        # 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 cost_raster is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_COST_RASTER))
        if cost_raster_band is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_RASTER_BAND))
        if polygons_to_connect is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_START_LAYER))
        if current_roads is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_START_LAYER))
        if skidding_distance is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_START_LAYER))
        if method_of_generation is None:
            raise QgsProcessingException(self.invalidSourceError(parameters, self.INPUT_START_LAYER))

        # We try to see if there are divergence between the CRSs of the inputs
        if cost_raster.crs() != polygons_to_connect.sourceCrs() \
                or polygons_to_connect.sourceCrs() != current_roads.sourceCrs():
            raise QgsProcessingException(self.tr("ERROR: The input layers have different CRSs."))

        # We check if the cost raster in indeed numeric
        if cost_raster.rasterType() not in [cost_raster.Multiband, cost_raster.GrayOrUndefined]:
            raise QgsProcessingException(self.tr("ERROR: The input cost raster is not numeric."))

        # We initialize the "sink", an object that will make use able to create an output.
        # First, we create the fields for the attributes of our lines as outputs.
        # They will only have one field :
        sink_fields = MinCostPathHelper.create_fields()
        # We indicate that our output will be a line, stored in WBK format.
        output_geometry_type = QgsWkbTypes.LineString
        # Finally, we create the field object and register the destination ID of it.
        # This sink will be based on the OUTPUT parameter we registered during initialization,
        # will have the fields and the geometry type we just created, and the same CRS as the cost raster.
        (sink, dest_id) = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            fields=sink_fields,
            geometryType=output_geometry_type,
            crs=cost_raster.crs(),
        )

        # If sink was not created, 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 invalidSinkError method to return a standard
        # helper text for when a sink cannot be evaluated
        if sink is None:
            raise QgsProcessingException(self.invalidSinkError(parameters, self.OUTPUT))

        # We check if the raster has been read correctly
        feedback.pushInfo(self.tr("The extent of the raster is : " + cost_raster.dataProvider().extent().asWktCoordinates()))
        # We put the data of the raster into a variable that we will send to the algorithm.
        block = MinCostPathHelper.get_all_block(cost_raster, cost_raster_band)
        # We transform the raster data into a matrix and check if the matrix contains negative values
        # CAREFUL : The matrix is created in a raster coordinate systems; rows (y axis) start at the top
        # and go to the bottom. This implies a transformation when getting back the values from a cartesian
        # system (rows go from bottom to top)
        matrix, contains_negative = MinCostPathHelper.block2matrix(block)
        # We display a feedback on the loading of the raster, or we display an error if needed
        if block.height() == 0 or block.width() == 0:
            raise QgsProcessingException(self.tr("ERROR: The raster couldn't be read properly (0 rows or 0 columns). "
                                                 "This is often due to the raster being too big. "
                                                 "Try to lower the resolution of your raster, and/or limit it to the"
                                                 "extent of your data."))
        feedback.pushInfo(self.tr("The size of the cost raster is: %d * %d pixels") % (block.height(), block.width()))

        # If there are negative values in the raster, we make an issue.
        if contains_negative:
            raise QgsProcessingException(self.tr("ERROR: Cost raster contains negative value."))

        feedback.pushInfo("Scanning the polygons to reach...")
        # First of all : We transform the starting polygons into cells on the raster (coordinates
        # in rows and colons).
        polygons_to_reach_features = list(polygons_to_connect.getFeatures())
        # feedback.pushInfo(str(len(start_features)))
        # We make a set of nodes to reach.
        set_of_nodes_to_reach, heuristicDictionnary = MinCostPathHelper.features_to_row_cols(polygons_to_reach_features,
                                                                                            heuristic_in_polygons_index,
                                                                                            cost_raster)
        # If there are no nodes to reach (e.g. all polygons are out of the raster)
        if len(set_of_nodes_to_reach) == 0:
            raise QgsProcessingException(self.tr("ERROR: There is no polygon to reach in this raster. Check if some"
                                                 "polygons are inside the raster."))
        feedback.pushInfo("Polygons scanned !")

        feedback.pushInfo("Scanning the existing roads...")
        # We do another set concerning the nodes that contains roads to connect to
        roads_to_connect_to_features = list(current_roads.getFeatures())
        # feedback.pushInfo(str(len(end_features)))
        set_of_nodes_to_connect_to, uselessHeuristicDictionary = MinCostPathHelper.features_to_row_cols(roads_to_connect_to_features,
                                                                                                        None,
                                                                                                        cost_raster)
        # If there is no nodes to connect to, throw an exception
        if len(set_of_nodes_to_connect_to) == 0:
            raise QgsProcessingException(self.tr("ERROR: There is no road to connect to in this raster. Check if some"
                                                 "roads are inside the raster."))
        # If some overlap, raise another exception :
        if set_of_nodes_to_reach in set_of_nodes_to_connect_to:
            raise QgsProcessingException(self.tr("ERROR: Some polygons to reach are overlapping with roads "
                                                 "to connect to given this resolution."))
        feedback.pushInfo("Roads scanned !")

        # Before we start, we need to order the nodes with the chosen heuristic.
        # We create a list that we are going to order.
        list_of_nodes_to_reach = list(set_of_nodes_to_reach)

        # If the method of generation asks for a random order, we shuffle the list randomly and it's over.
        if method_of_generation == '0':
            feedback.pushInfo("Randomizing order of cells to visit...")
            random.shuffle(list_of_nodes_to_reach)
        # If not, we create a list that will contain the minimal distance between the given node and the nodes to
        # connect to.
        else:
            list_of_nodes_to_reach_with_order = list()

            feedback.pushInfo("Computing distances between polygons and roads...(This can take some time !)")
            feedbackProgress = 0

            # To quickly calculate the distance from the existing roads to each node in our polygons, we will use
            # the k-d tree function from Scipy.
            # For that, we need to make an Numpy array containing our road nodes
            pointsToReach = list()
            for node in set_of_nodes_to_connect_to:
                nodeToPoint = MinCostPathHelper._row_col_to_point(node, cost_raster)
                pointsToReach.append(nodeToPoint)
            numpyArrayOfRoadPoints = np.array(pointsToReach)
            spatialKDTREEForDistanceSearch = KDTree(numpyArrayOfRoadPoints, leafsize=20)

            # Then, we search and find the distance from current roads for each point.
            for node in list_of_nodes_to_reach:
                # feedback.pushInfo("For node " + str(node) + "...")

                # Obsolete with the use of the k-d tree.
                # minimalDistance = MinCostPathHelper.minimum_distance_to_a_node(node,
                                                                               # set_of_nodes_to_connect_to,
                                                                               # cost_raster)
                # feedback.pushInfo("Brute force gives minimal distance : " + str(minimalDistance))

                # We now make a query to the tree to find the closest node of the tree to the given node.
                nodeAsPoint = MinCostPathHelper._row_col_to_point(node, cost_raster)
                minimalDistance = spatialKDTREEForDistanceSearch.query(nodeAsPoint)[0]
                # feedback.pushInfo("K-D Tree gives minimal distance : " + str(minimalDistance))

                # We add the distance to a tuple, so that we can classify the node afterward.
                list_of_nodes_to_reach_with_order.append((minimalDistance, node))
                feedbackProgress += 1
                feedback.setProgress(100 * (feedbackProgress / len(list_of_nodes_to_reach)))
                if feedback.isCanceled():
                    raise QgsProcessingException(self.tr("ERROR: Operation was cancelled."))

            feedback.pushInfo("Computing distances is done !")

            # We then sort according to this distance, in increasing or decreasing order.
            if method_of_generation == '1':
                # We sort the list by the heuristic inside the polygons, AND the main heuristic.
                list_of_nodes_to_reach_with_order = sorted(list_of_nodes_to_reach_with_order,
                                                           key=lambda tuple: (heuristicDictionnary[tuple[1]],
                                                                              tuple[0]))
                # list_of_nodes_to_reach_with_order.sort()
                feedback.pushInfo("Ordering towards closest cells to visit...")
            else:
                feedback.pushInfo("Ordering towards farthest cells to visit...")
                # Here, we take the opposite of the distance so that it respects the "smallest first" sorting, but
                # in the opposite way to select the farthest first.
                list_of_nodes_to_reach_with_order = sorted(list_of_nodes_to_reach_with_order,
                                                           key=lambda tuple: (heuristicDictionnary[tuple[1]],
                                                                              -tuple[0]))
                # list_of_nodes_to_reach_with_order.sort(reverse=True)

            feedback.pushInfo("Ordering is done !")

            # We put the result in the list of nodes to reach back again, removing the distance.
            list_of_nodes_to_reach = [i[1] for i in list_of_nodes_to_reach_with_order]

        # Now, time to launch the algorithm properly !
        feedback.pushInfo(self.tr("Generating the road network...(This can take some time !)"))

        # First, we have to initialize some things.
        feedbackProgress = 0
        errorMessages = 0
        listOfResults = list()
        pointsToReach = set()
        for node in set_of_nodes_to_connect_to:
            nodeToPoint = MinCostPathHelper._row_col_to_point(node, cost_raster)
            pointsToReach.add(nodeToPoint)

        # This is the road matrix. It's use to quickly know if there are roads at two given coordinates.
        roadMatrix = np.zeros( (cost_raster.height(), cost_raster.width()) )
        # A "1" means that there is a road at the given coordinate.
        for node in set_of_nodes_to_connect_to:
            roadMatrix[node[0]][node[1]] = 1

        # Now, we have to initialize the circle neighborhood.
        skiddingDistanceCircleNeighborhood = MinCostPathHelper.createRelativeCircleNeighborhood(skidding_distance,
                                                                                                cost_raster)
        feedback.pushInfo(self.tr("Size of the skidding neighborhood : " + str(len(skiddingDistanceCircleNeighborhood))))

        for nodeToReach in list_of_nodes_to_reach:
            feedbackProgress += 1

            # If the node to reach is inside a no-value pixel, no need to look at it.
            if matrix[(len(matrix)-1)-nodeToReach[0]][nodeToReach[1]] is not None:
                # First, we check the distance between the node and the nodes to connect to,
                # to see if it's not at a skidding distance of it.

                # Obsolete with the k-d tree.
                # minimalDistanceToNodesToConnect = MinCostPathHelper.minimum_distance_to_a_node(nodeToReach,
                                                                                               # set_of_nodes_to_connect_to,
                                                                                               # cost_raster)

                # To quickly calculate the distance from the existing roads to each node in our polygons, we will use
                # the k-d tree function from SciPy.
                # Now obsolete with the circle neighborhood method.
                # For that, we need to make an Numpy array containing our road nodes
                # spatialKDTREEForDistanceSearch = KDTree(np.array(list(pointsToReach)), leafsize=20)
                # nodeAsPoint = MinCostPathHelper._row_col_to_point(nodeToReach, cost_raster)
                # minimalDistanceToNodesToConnect = spatialKDTREEForDistanceSearch.query(nodeAsPoint)[0]

                # If it's superior, we create a road to this node
                # if minimalDistanceToNodesToConnect > skidding_distance:

                # New method : using a relative neighborhood.
                foundRoadInSkiddingDistance = MinCostPathHelper.checkRelativeCircleNeighborhoodForRoads(skiddingDistanceCircleNeighborhood,
                                                                                                        nodeToReach,
                                                                                                        roadMatrix)

                if not foundRoadInSkiddingDistance :
                    start_row_col = nodeToReach
                    end_row_cols = list(set_of_nodes_to_connect_to)
                    min_cost_path, costs, selected_end = dijkstra(start_row_col, end_row_cols, matrix,
                                                                  angles_considered, punisherAngleDictionnary, feedback)
                    # If there was a problem, we indicate if it's because the search was cancelled by the user
                    # or if there was no end point that could be reached.
                    if min_cost_path is None:
                        if feedback.isCanceled():
                            raise QgsProcessingException(self.tr("ERROR: Search canceled."))
                        else:
                            errorMessages += 1

                    # If there wasn't a problem, we save the results
                    else:
                        # When the road is done by the Dijkstra algorithm, we put the path and the cost
                        # in the list of results
                        listOfResults.append((min_cost_path, costs[-1]))
                        # We also add the nodes of the created path to the set of nodes that can be reached now
                        set_of_nodes_to_connect_to.update(min_cost_path)
                        for node in min_cost_path:
                            nodeToPoint = MinCostPathHelper._row_col_to_point(node, cost_raster)
                            pointsToReach.add(nodeToPoint)
                            roadMatrix[node[0]][node[1]] = 1

            feedback.setProgress(100 * (feedbackProgress / len(list_of_nodes_to_reach)))

        # When the loop is done..
        feedback.setProgress(100)
        feedback.pushInfo(self.tr("Network created ! Saving network..."))

        # For every path we create, we save it as a line and put it into the sink !
        ID = 1
        for (path, cost) in listOfResults:
            # feedback.pushInfo("Cost of feature saved : " + str(cost))
            # Time to save the path as a vector.
            # We take the starting and ending points as pointXY
            start_point = MinCostPathHelper._row_col_to_point(path[0], cost_raster)
            end_point = MinCostPathHelper._row_col_to_point(path[-1], cost_raster)
            # We make a list of Qgs.pointXY from the nodes in our pathlist
            path_points = MinCostPathHelper.create_points_from_path(cost_raster, path, start_point, end_point)
            # With the total cost which is the last item in our accumulated cost list,
            # we create the PolyLine that will be returned as a vector.
            path_feature = MinCostPathHelper.create_path_feature_from_points(path_points, cost, ID, sink_fields)
            # Into the sink that serves as our output, we put the PolyLines from the list of lines we created
            # one by one
            sink.addFeature(path_feature, QgsFeatureSink.FastInsert)
            ID += 1

        # We display the error messages if there was some
        if errorMessages > 0:
            feedback.pushInfo("WARNING : During the pathfinding, there was " + str(errorMessages) + " cases were a road could not"
                              " be constructed to a certain point to reach in a harvested polygon. When this happens, it's often due to"
                              " the cost raster containing No Data pixels that the algorithm cannot cross. Please, check"
                              " your cost raster and and change the No Data pixels if needed.")

        # When all is done, we return our output that is linked to the sink.
        return {self.OUTPUT: dest_id}

    # 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 'Forest Road Network Creation'

    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 ForestRoadNetworkAlgorithm()

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

    def shortHelpString(self):
        return self.tr("""
        This algorithm creates a forest road network based on areas to access (polygons) and current roads to connect them to (lines).
        
        **Parameters:**
          
          Please ensure all the input layers have the same CRS.
        
          - Cost raster layer: Numeric raster layer that represents the cost of each spatial unit. It should not contains negative value. Pixel with `NoData` value represent it is unreachable.
         
          - Cost raster band: The input band of the cost raster.
         
          - Polygons to access via the generated roads: Layer that contains the polygons to access.
         
          - Roads to connect the polygons to access to: Layer that contains the roads to connect the polygons to.
          
          - Skidding distance. Maximum distance that a cell can be to not need a road going up to it.
          
          - Method of generation : a parameter indicating what type of heuristic is used to generate the network. Random cell order, farther cells from current roads first, closer cells from curent roads first.
          
          - Attribute containing an heuristic : An attribute field of the polygons that contains an heuristic that describe in which order the algorithm should reach them. The lower the value, the higher the priority; this way, the heuristic can be a date or a time. It is combined with the heuristic chosen before by the user to determine the order in which pixels are accessed a single polygon.
         
          - Punishment of angles : The algorithm can "punish" the use of steep angle between a pixel and another by increasing the cost when it is used to allow for smoother lines. Remember that this does not affect the angle of connection between the roads, but only the zigzags in each road.
         
        """)

    def shortDescription(self):
        return self.tr('Generate a network of roads to connect forest area to an existing road network.')

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

    def tags(self):
        return ['least', 'cost', 'path', 'distance', 'raster', 'analysis', 'road', 'network', 'forest', 'A*', 'dijkstra']

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

    # Function to transform a given row/column into a QGIS point with a x,y
    # coordinates based on the resolution of the raster layer we're considering
    # (calculated with its extent and number of cells)
    # CAREFUL : Coordinate system is cartesian, in opposition to the
    # matrix containing the data of the raster (see further down)
    @staticmethod
    def _row_col_to_point(row_col, raster_layer):
        xres = raster_layer.rasterUnitsPerPixelX()
        yres = raster_layer.rasterUnitsPerPixelY()
        extent = raster_layer.dataProvider().extent()

        x = (row_col[1] + 0.5) * xres + extent.xMinimum()
        # There is a dissonance about how I see y axis of the raster
        # and how the dijstra algorithm sees it. This is where the transformation is made.
        y = (row_col[0] + 0.5) * yres + extent.yMinimum()
        return QgsPointXY(x, y)

    # Method to determine where a given polygon is in the raster
    @staticmethod
    def _polygon_to_row_col(polygon, raster_layer):
        # We initialize the object to return : a set of nodes, corresponding a tuple of a row
        # and a column in the raster
        listOfNodesInPolygons = set()

        # We get the extent of the raster
        xres = raster_layer.rasterUnitsPerPixelX()
        yres = raster_layer.rasterUnitsPerPixelY()
        extentRaster = raster_layer.dataProvider().extent()
        maxRasterCols = round((extentRaster.xMaximum() - extentRaster.xMinimum()) / xres)
        maxRasterRows = round((extentRaster.yMaximum() - extentRaster.yMinimum()) / yres)

        # We get the extent of the polygon
        extentPolygon = QgsGeometry.fromPolygonXY(polygon).boundingBox()

        # We traduce this extent into limits of columns and rows on the raster
        rowMin = floor((extentPolygon.yMinimum() - extentRaster.yMinimum()) / yres)
        rowMax = floor((extentPolygon.yMaximum() - extentRaster.yMinimum()) / yres)
        colMin = floor((extentPolygon.xMinimum() - extentRaster.xMinimum()) / xres)
        colMax = floor((extentPolygon.xMaximum() - extentRaster.xMinimum()) / xres)

        # If one of these values is not in the range of the raster, then we
        # restrict it to column and rows that are inside of it.
        if rowMin < 0: rowMin = 0
        if rowMax > maxRasterRows: rowMax = maxRasterRows
        if colMin < 0: colMin = 0
        if colMax > maxRasterCols: colMax = maxRasterCols
        # If the polygon is out of range, then we return an empty set
        if rowMin > maxRasterRows or colMin > maxRasterCols or rowMax < 0 or colMax < 0:
            return set()

        # If not, we treat the selected range of the raster to avoid looking
        # at all of the raster.
        # For each column of our range,
        for col in range(colMin, colMax):
            # For each row of our range,
            for row in range(rowMin, rowMax):
                # We make a centroid of the row/column of the raster
                centroid = MinCostPathHelper._row_col_to_point((row, col), raster_layer)
                # We check if the centroid is in the polygon
                isInPolygon = QgsGeometry.fromPolygonXY(polygon).contains(centroid)

                # If it is, we add the centroid to a set with the form of a node (tuple (row, column))
                if isInPolygon:
                    listOfNodesInPolygons.add((row, col))
                # if not, we continue the loops

        # When the loops are done, we return the set of nodes that have been taken into account
        return listOfNodesInPolygons

    # Method to determine where a given line is in the raster. Similar to previous.
    @staticmethod
    def _line_to_row_col(line, raster_layer):
        # We initialize the object to return : a set of nodes, corresponding a tuple of a row
        # and a column in the raster
        listOfNodesInPolygons = set()

        # We get the extent of the raster
        xres = raster_layer.rasterUnitsPerPixelX()
        yres = raster_layer.rasterUnitsPerPixelY()
        extentRaster = raster_layer.dataProvider().extent()
        maxRasterCols = round((extentRaster.xMaximum() - extentRaster.xMinimum()) / xres)
        maxRasterRows = round((extentRaster.yMaximum() - extentRaster.yMinimum()) / yres)

        # We get the extent of the line
        extentLine = QgsGeometry.fromPolylineXY(line).boundingBox()

        # We traduce this extent into limits of colomns and rows on the raster
        rowMin = floor((extentLine.yMinimum() - extentRaster.yMinimum()) / yres)
        rowMax = floor((extentLine.yMaximum() - extentRaster.yMinimum()) / yres)
        colMin = floor((extentLine.xMinimum() - extentRaster.xMinimum()) / xres)
        colMax = floor((extentLine.xMaximum() - extentRaster.xMinimum()) / xres)

        # If one of these values is not in the range of the raster, then we
        # restrict it to column and rows that are inside of it.
        if rowMin < 0: rowMin = 0
        if rowMax > maxRasterRows: rowMax = maxRasterRows
        if colMin < 0: colMin = 0
        if colMax > maxRasterCols: colMax = maxRasterCols
        # If the polygon is out of range, then we return an empty set
        if rowMin > maxRasterRows or colMin > maxRasterCols or rowMax < 0 or colMax < 0:
            return set()

        # If not, we treat the select range of the raster.
        # For each column of our range,
        for col in range(colMin, colMax):
            # For each row of our range,
            for row in range(rowMin, rowMax):
                # We make square polygon around the cell
                centroid = MinCostPathHelper._row_col_to_point((row, col), raster_layer)
                halfACellX = 0.5 * raster_layer.rasterUnitsPerPixelX()
                halfACellY = 0.5 * raster_layer.rasterUnitsPerPixelY()
                square = QgsGeometry.fromPolygonXY([[QgsPointXY(centroid.x() - halfACellX, centroid.y() - halfACellY),
                                                   QgsPointXY(centroid.x() + halfACellX, centroid.y() - halfACellY),
                                                   QgsPointXY(centroid.x() + halfACellX, centroid.y() + halfACellY),
                                                   QgsPointXY(centroid.x() - halfACellX, centroid.y() + halfACellY)]])

                # We check if the centroid is in the polygon
                intersectWithSquare = square.intersects(QgsGeometry.fromPolylineXY(line))

                # If it is, we add the centroid to a set with the form of a node (tuple (row, column))
                if intersectWithSquare:
                    listOfNodesInPolygons.add((row, col))
                # if not, we continue the loops

        # When the loops are done, we return the set of nodes that have been taken into account
        return listOfNodesInPolygons

    # Function to return a list of Qgs.pointXY. Each point is made based on the center of the node
    # that we get from the path list.
    # At the end, we put the precise coordinates of the starting/ending nodes that were given by
    # the user at the start.
    @staticmethod
    def create_points_from_path(cost_raster, min_cost_path, start_point, end_point):
        path_points = list(
            map(lambda row_col: MinCostPathHelper._row_col_to_point(row_col, cost_raster), min_cost_path))
        path_points[0].setX(start_point.x())
        path_points[0].setY(start_point.y())
        path_points[-1].setX(end_point.x())
        path_points[-1].setY(end_point.y())
        return path_points

    @staticmethod
    def create_fields():
        # Create an ID field to know in which order the roads have been constructed
        id_field = QgsField("Construction order", QVariant.Int, "integer", 10, 3)
        # Create the field of "total cost" by indicating name, type, typeName,
        # lenght and precision (decimals in that case)
        cost_field = QgsField("Total cost", QVariant.Double, "double", 15, 3)
        # Then, we create a container of multiple fields
        fields = QgsFields()
        # We add the fields to the container
        fields.append(id_field)
        fields.append(cost_field)
        # We return the container with our fields.
        return fields

    # Function to create a polyline with the list of qgs.pointXY
    @staticmethod
    def create_path_feature_from_points(path_points, total_cost, ID, fields):
        # We create the geometry of the polyline
        polyline = QgsGeometry.fromPolylineXY(path_points)
        # We retrieve the fields and add them to the feature
        feature = QgsFeature(fields)
        cost_index = feature.fieldNameIndex("total cost")
        feature.setAttribute(cost_index, total_cost)  # cost
        id_index = feature.fieldNameIndex("Construction order")
        feature.setAttribute(id_index, ID) # id
        # We add the geometry to the feature
        feature.setGeometry(polyline)
        return feature

    # Method to transform given features into a set of
    # nodes (row + column) on the raster.
    # Features have to be lines or polygons.
    # Also return a dictionary containing the heuristic read in the polygon or line
    # for use in ordering the nodes to reach.
    @staticmethod
    def features_to_row_cols(given_features, heuristic_index, raster_layer):

        row_cols = set()
        heuristic_dictionary = dict()
        # extent = raster_layer.dataProvider().extent()
        # if extent.isNull() or extent.isEmpty:
        #     return list(col_rows)

        for given_feature in given_features:
            if given_feature.hasGeometry():
                given_feature_geom = given_feature.geometry()

                if heuristic_index is not None:
                    attributes = given_feature.attributes()
                    heuristic = attributes[heuristic_index]
                else:
                    heuristic = 0

                # Case of multipolygons
                if given_feature_geom.wkbType() == QgsWkbTypes.MultiPolygon:
                    multi_polygon = given_feature_geom.asMultiPolygon()
                    for polygon in multi_polygon:
                        row_cols_for_this_polygon = MinCostPathHelper._polygon_to_row_col(polygon, raster_layer)
                        row_cols.update(row_cols_for_this_polygon)
                        for row_col in row_cols_for_this_polygon:
                            heuristic_dictionary[row_col] = heuristic

                # Case of polygons
                elif given_feature_geom.wkbType() == QgsWkbTypes.Polygon:
                    polygon = given_feature_geom.asPolygon()
                    row_cols_for_this_polygon = MinCostPathHelper._polygon_to_row_col(polygon, raster_layer)
                    row_cols.update(row_cols_for_this_polygon)
                    for row_col in row_cols_for_this_polygon:
                        heuristic_dictionary[row_col] = heuristic

                # Case of multi lines
                elif given_feature_geom.wkbType() == QgsWkbTypes.MultiLineString:
                    multi_line = given_feature_geom.asMultiPolyline()
                    for line in multi_line:
                        row_cols_for_this_line = MinCostPathHelper._line_to_row_col(line, raster_layer)
                        row_cols.update(row_cols_for_this_line)
                        for row_col in row_cols_for_this_line:
                            heuristic_dictionary[row_col] = heuristic

                # Case of lines
                elif given_feature_geom.wkbType() == QgsWkbTypes.LineString:
                    line = given_feature_geom.asPolyline()
                    row_cols_for_this_line = MinCostPathHelper._line_to_row_col(line, raster_layer)
                    row_cols.update(row_cols_for_this_line)
                    for row_col in row_cols_for_this_line:
                        heuristic_dictionary[row_col] = heuristic

        return row_cols, heuristic_dictionary

    # Function that get the data block from a entire raster for a given band
    @staticmethod
    def get_all_block(raster_layer, band_num):
        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(band_num, extent, width, height)

    # Function that transforms
    @staticmethod
    def block2matrix(block):
        contains_negative = False
        matrix = [[None if block.isNoData(i, j) else block.value(i, j) for j in range(block.width())]
                  for i in range(block.height())]

        for l in matrix:
            for v in l:
                if v is not None:
                    if v < 0:
                        contains_negative = True

        return matrix, contains_negative

    # This function return the minimum distance between a given node, and the nodes in a set or list of nodes.
    @staticmethod
    def minimum_distance_to_a_node(node, listOrSetOfNodes, raster_layer):

        listOfNodes = list(listOrSetOfNodes)
        minimumDistance = float("inf")
        pointGeometryOfNode = QgsGeometry.fromPointXY(MinCostPathHelper._row_col_to_point(node, raster_layer))

        for otherNode in listOfNodes:
            distance = pointGeometryOfNode.distance(QgsGeometry.fromPointXY(MinCostPathHelper._row_col_to_point(otherNode, raster_layer)))

            if distance < minimumDistance:
                minimumDistance = distance

        return minimumDistance

    @staticmethod
    def createRelativeCircleNeighborhood(skiddingDistance, raster_layer):
        """"This method initialize a relative circle neighborhood based on the size of the pixels and on the
        skidding distance inputted by the user, in order to check rapidly if an existing road is at skidding distance
        from a given node."""

        xres = raster_layer.rasterUnitsPerPixelX()
        yres = raster_layer.rasterUnitsPerPixelY()

        widthOfNeighborhood = floor(skiddingDistance / xres) + 1
        heightOfNeighborhood = floor(skiddingDistance / yres) + 1

        relativeCircleNeighborhood = set()

        for col in range(-widthOfNeighborhood, widthOfNeighborhood+1):
            for row in range(-heightOfNeighborhood, heightOfNeighborhood+1):
                # If the euclidian distance between a relative coordinate and the "origin" cell is superior to the
                # skidding distance, this cell won't be part of the neighborhood.
                if sqrt((row*yres)**2 + (col*xres)**2) <= skiddingDistance:
                    relativeCircleNeighborhood.add((row, col))

        return relativeCircleNeighborhood

    @staticmethod
    def checkRelativeCircleNeighborhoodForRoads(relativeCircleNeighborhood, nodeAsRowCol, roadMatrix):
        """"This function uses the relative circle neighborhood to check if a road is at skidding distance from a
        node."""

        rowOfNode = nodeAsRowCol[0]
        colOfNode = nodeAsRowCol[1]
        foundARoad = False

        for relativeNeighbour in relativeCircleNeighborhood:
            relativeRow = rowOfNode + relativeNeighbour[0]
            relativeColumn = colOfNode + relativeNeighbour[1]
            # We check if the relative row index is correct; same for the column.
            if relativeRow > 0 and relativeRow < len(roadMatrix):
                if relativeColumn > 0 and relativeColumn < len(roadMatrix[0]):
                    if roadMatrix[relativeRow][relativeColumn] == 1:
                        foundARoad = True
                        break

        return foundARoad