# -*- coding: utf-8 -*-
"""
/***************************************************************************
 DsgTools
                                 A QGIS plugin
 Brazilian Army Cartographic Production Tools
                              -------------------
        begin                : 2023-05-19
        git sha              : $Format:%H$
        copyright            : (C) 2023 by Philipe Borba - Cartographic Engineer @ Brazilian Army
        email                : borba.philipe@eb.mil.br
 ***************************************************************************/
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
"""

from collections import defaultdict
import math
from typing import Dict, List, Optional, Tuple, Union
from uuid import uuid4
import numpy as np
from itertools import islice
from DsgTools.core.DSGToolsProcessingAlgs.algRunner import AlgRunner
from DsgTools.core.GeometricTools import geometryHandler, rasterHandler
from DsgTools.core.GeometricTools.affine import Affine
from DsgTools.core.GeometricTools.layerHandler import LayerHandler
from DsgTools.core.GeometricTools.spatialRelationsHandler import SpatialRelationsHandler
from qgis.PyQt.QtCore import QCoreApplication, QMetaType, QByteArray
from qgis.core import (
    QgsFeature,
    QgsFeatureSink,
    QgsField,
    QgsFields,
    QgsProcessing,
    QgsProcessingAlgorithm,
    QgsProcessingMultiStepFeedback,
    QgsProcessingParameterFeatureSink,
    QgsProcessingParameterNumber,
    QgsProcessingParameterRasterLayer,
    QgsProcessingUtils,
    QgsVectorLayer,
    QgsWkbTypes,
    QgsProcessingParameterEnum,
    QgsCoordinateReferenceSystem,
    QgsRasterLayer,
    QgsProcessingContext,
    QgsFeedback,
    QgsProcessingParameterField,
    QgsVectorLayerUtils,
    QgsProcessingParameterVectorLayer,
    QgsProcessingParameterExpression,
    QgsProcessingParameterBoolean,
    QgsFeatureRequest,
    QgsProcessingParameterCrs,
    Qgis,
    QgsProcessingParameterDefinition,
    QgsCoordinateTransform,
    QgsProject,
)


class ExtractElevationPoints(QgsProcessingAlgorithm):

    INPUT_DEM = "INPUT_DEM"
    CONTOUR_LINES = "CONTOUR_LINES"
    CONTOUR_ATTR = "CONTOUR_ATTR"
    DEPRESSION_EXPRESSION = "DEPRESSION_EXPRESSION"
    CONTOUR_INTERVAL = "CONTOUR_INTERVAL"
    GEOGRAPHIC_BOUNDARY = "GEOGRAPHIC_BOUNDARY"
    SCALE = "SCALE"
    WATER_BODIES = "WATER_BODIES"
    AREA_WITHOUT_INFORMATION_POLYGONS = "AREA_WITHOUT_INFORMATION_POLYGONS"
    NATURAL_POINT_FEATURES = "NATURAL_POINT_FEATURES"
    DRAINAGE_LINES = "DRAINAGE_LINES"
    DRAINAGE_LINES_NAME_ATTRIBUTE = "DRAINAGE_LINES_NAME_ATTRIBUTE"
    CONSIDER_DRAINAGES_WITHOUT_NAMES = "CONSIDER_DRAINAGES_WITHOUT_NAMES"
    ROADS = "ROADS"
    MAIN_ROADS_EXPRESSION = "MAIN_ROADS_EXPRESSION"
    OTHER_ROADS_EXPRESSION = "OTHER_ROADS_EXPRESSION"
    ONLY_HILLTOPS = "ONLY_HILLTOPS"
    OUTPUT = "OUTPUT"

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INPUT_DEM,
                self.tr("Input DEM"),
            )
        )
        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.CONTOUR_LINES,
                self.tr("Contour Lines"),
                [QgsProcessing.TypeVectorLine],
                defaultValue="elemnat_curva_nivel_l",
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.CONTOUR_ATTR,
                self.tr("Contour value field"),
                "cota",
                self.CONTOUR_LINES,
                QgsProcessingParameterField.Numeric,
            )
        )

        self.addParameter(
            QgsProcessingParameterExpression(
                self.DEPRESSION_EXPRESSION,
                self.tr("Filter expression for contour that are depressions."),
                """ "depressao" = 1 """,
                self.CONTOUR_LINES,
                optional=True,
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.CONTOUR_INTERVAL,
                self.tr("Contour interval"),
                minValue=0,
                defaultValue=10,
            )
        )

        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.GEOGRAPHIC_BOUNDARY,
                self.tr("Geographic bounds layer"),
                [QgsProcessing.TypeVectorPolygon],
                optional=False,
                defaultValue="aux_moldura_a",
            )
        )

        self.scales = [
            "1:25.000",
            "1:50.000",
            "1:100.000",
            "1:250.000",
        ]

        self.distances = {
            0: 20e-3 * 25_000,
            1: 20e-3 * 50_000,
            2: 20e-3 * 100_000,
            3: 20e-3 * 250_000,
        }

        self.minContourLenghts = {
            0: 8e-3 * 25_000,
            1: 8e-3 * 50_000,
            2: 8e-3 * 100_000,
            3: 8e-3 * 250_000,
        }
        self.gridSpacingDict = {
            0: 5000,
            1: 10000,
            2: 20000,
            3: 40000,
        }
        self.contourBufferLengths = {
            0: 6.4e-9 * 25_000,
            1: 6.4e-9 * 50_000,
            2: 6.4e-9 * 100_000,
            3: 6.4e-9 * 250_000,
        }

        self.planeGridSpacingDict = {
            0: 20e-3 * 25_000,
            1: 20e-3 * 50_000,
            2: 20e-3 * 100_000,
            3: 20e-3 * 250_000,
        }

        self.addParameter(
            QgsProcessingParameterEnum(
                self.SCALE, self.tr("Scale"), options=self.scales, defaultValue=0
            )
        )

        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.WATER_BODIES,
                self.tr("Water Bodies"),
                [QgsProcessing.TypeVectorPolygon],
                defaultValue="cobter_massa_dagua_a",
            )
        )

        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.AREA_WITHOUT_INFORMATION_POLYGONS,
                self.tr("Area without information layer"),
                [QgsProcessing.TypeVectorPolygon],
                optional=True,
            )
        )

        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.DRAINAGE_LINES,
                self.tr("Drainage lines"),
                [QgsProcessing.TypeVectorLine],
                defaultValue="elemnat_trecho_drenagem_l",
            )
        )

        self.addParameter(
            QgsProcessingParameterField(
                self.DRAINAGE_LINES_NAME_ATTRIBUTE,
                self.tr("Drainage name field"),
                "nome",
                self.DRAINAGE_LINES,
                QgsProcessingParameterField.String,
            )
        )

        self.addParameter(
            QgsProcessingParameterBoolean(
                self.CONSIDER_DRAINAGES_WITHOUT_NAMES,
                self.tr("Consider drainages without names"),
                defaultValue=True,
            )
        )

        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.ROADS,
                self.tr("Roads"),
                [QgsProcessing.TypeVectorLine],
                defaultValue="infra_via_deslocamento_l",
            )
        )

        self.addParameter(
            QgsProcessingParameterExpression(
                self.MAIN_ROADS_EXPRESSION,
                self.tr("Filter expression for main roads"),
                """"tipo" = 4 or "jurisdicao" in (1,2)""",
                self.ROADS,
            )
        )

        self.addParameter(
            QgsProcessingParameterExpression(
                self.OTHER_ROADS_EXPRESSION,
                self.tr("Filter expression for other roads"),
                """"tipo" in (2,4)""",
                self.ROADS,
                optional=True,
            )
        )

        self.addParameter(
            QgsProcessingParameterBoolean(
                self.ONLY_HILLTOPS,
                self.tr(
                    "Extract only hilltops (extracts all hilltop points, "
                    "including depressions, and discards all other methods)"
                ),
                defaultValue=False,
            )
        )

        self.addParameter(
            QgsProcessingParameterFeatureSink(self.OUTPUT, "Output spot elevation")
        )

    def processAlgorithm(self, parameters, context, feedback):
        self.algRunner = AlgRunner()
        inputRaster = self.parameterAsRasterLayer(parameters, self.INPUT_DEM, context)
        contourLyr = self.parameterAsVectorLayer(
            parameters, self.CONTOUR_LINES, context
        )
        heightFieldName = self.parameterAsStrings(
            parameters, self.CONTOUR_ATTR, context
        )[0]
        contourHeightInterval = self.parameterAsDouble(
            parameters, self.CONTOUR_INTERVAL, context
        )
        depressionExpression = self.parameterAsExpression(
            parameters, self.DEPRESSION_EXPRESSION, context
        )
        if depressionExpression == "":
            depressionExpression = None
        geographicBoundaryLyr = self.parameterAsVectorLayer(
            parameters, self.GEOGRAPHIC_BOUNDARY, context
        )
        scale = self.parameterAsEnum(parameters, self.SCALE, context)
        areaWithoutInformationLyr = self.parameterAsVectorLayer(
            parameters, self.AREA_WITHOUT_INFORMATION_POLYGONS, context
        )
        waterBodiesLyr = self.parameterAsVectorLayer(
            parameters, self.WATER_BODIES, context
        )
        naturalPointFeaturesLyr = self.parameterAsVectorLayer(
            parameters, self.NATURAL_POINT_FEATURES, context
        )
        drainagesLyr = self.parameterAsVectorLayer(
            parameters, self.DRAINAGE_LINES, context
        )
        drainageNameAttribute = self.parameterAsStrings(
            parameters, self.DRAINAGE_LINES_NAME_ATTRIBUTE, context
        )[0]
        getDrainagesWithoutName = self.parameterAsBool(
            parameters, self.CONSIDER_DRAINAGES_WITHOUT_NAMES, context
        )

        roadsLyr = self.parameterAsVectorLayer(parameters, self.ROADS, context)
        mainRoadsFilterExpression = self.parameterAsExpression(
            parameters, self.MAIN_ROADS_EXPRESSION, context
        )
        if mainRoadsFilterExpression == "":
            mainRoadsFilterExpression = None

        otherRoadsFilterExpression = self.parameterAsExpression(
            parameters, self.OTHER_ROADS_EXPRESSION, context
        )
        if otherRoadsFilterExpression == "":
            otherRoadsFilterExpression = None

        onlyHilltops = self.parameterAsBool(parameters, self.ONLY_HILLTOPS, context)

        self.outputCrs = QgsCoordinateReferenceSystem("EPSG:4674")

        self.bufferDist = self.distances[scale]
        self.minContourLength = self.minContourLenghts[scale]
        self.gridSpacing = self.gridSpacingDict[scale]
        self.planeGridSpacing = self.gridSpacingDict[scale]
        self.contourBufferLength = self.contourBufferLengths[scale]

        fields = QgsFields()
        fields.append(QgsField("cota", QMetaType.Type.Int))
        fields.append(QgsField("cota_mais_alta", QMetaType.Type.Int))
        fields.append(QgsField("cota_comprovada", QMetaType.Type.Int))
        fields.append(QgsField("ancora_horizontal", QMetaType.Type.Int))
        fields.append(QgsField("ancora_vertical", QMetaType.Type.Int))
        fields.append(QgsField("suprimir_simbologia", QMetaType.Type.Int))
        fields.append(QgsField("visivel", QMetaType.Type.Int))
        self.defaultAttrMap = {
            "cota_mais_alta": 2,
            "cota_comprovada": 2,
            "ancora_horizontal": 1,
            "ancora_vertical": 1,
            "suprimir_simbologia": 2,
            "visivel": 1,
        }

        (self.sink, self.sink_id) = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            fields,
            QgsWkbTypes.Point,
            self.outputCrs,
        )
        layerHandler = LayerHandler()
        nFeats = geographicBoundaryLyr.featureCount()
        nSteps = (
            nFeats
            + 7
            + 3 * (otherRoadsFilterExpression is not None)
            + 3 * getDrainagesWithoutName
        )
        multiStepFeedback = QgsProcessingMultiStepFeedback(nSteps, feedback)
        currentStep = 0
        if geographicBoundaryLyr.crs() != QgsCoordinateReferenceSystem("EPSG:4674"):
            geographicBoundaryLyr = self.algRunner.runReprojectLayer(
                layer=geographicBoundaryLyr,
                targetCrs=QgsCoordinateReferenceSystem("EPSG:4674"),
                feedback=multiStepFeedback,
                context=context,
            )
        currentStep += 1
        multiStepFeedback.setCurrentStep(currentStep)
        if inputRaster.crs() != QgsCoordinateReferenceSystem("EPSG:4674"):
            inputRaster = self.algRunner.runGdalWarp(
                rasterLayer=inputRaster,
                targetCrs=QgsCoordinateReferenceSystem("EPSG:4674"),
                context=context,
                feedback=multiStepFeedback,
            )

        currentStep += 1
        multiStepFeedback.setCurrentStep(currentStep)
        multiStepFeedback.pushInfo(self.tr("Getting roads"))
        mainRoadsLyr = (
            self.algRunner.runFilterExpression(
                inputLyr=roadsLyr,
                expression=mainRoadsFilterExpression,
                context=context,
                feedback=multiStepFeedback,
                is_child_algorithm=False,
            )
            if mainRoadsFilterExpression is not None
            else roadsLyr
        )
        currentStep += 1
        multiStepFeedback.setCurrentStep(currentStep)
        if mainRoadsLyr.crs() != QgsCoordinateReferenceSystem("EPSG:4674"):
            mainRoadsLyr = self.algRunner.runReprojectLayer(
                layer=mainRoadsLyr,
                targetCrs=QgsCoordinateReferenceSystem("EPSG:4674"),
                feedback=multiStepFeedback,
                context=context,
            )

        currentStep += 1
        multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            mainRoadsLyr, context, multiStepFeedback, is_child_algorithm=True
        )
        currentStep += 1

        if otherRoadsFilterExpression is not None:
            multiStepFeedback.setCurrentStep(currentStep)
            otherRoadsLyr = self.algRunner.runFilterExpression(
                inputLyr=roadsLyr,
                expression=otherRoadsFilterExpression,
                context=context,
                feedback=multiStepFeedback,
                is_child_algorithm=True,
            )
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            otherRoadsLyr = self.algRunner.runReprojectLayer(
                layer=otherRoadsLyr,
                targetCrs=QgsCoordinateReferenceSystem("EPSG:4674"),
                feedback=multiStepFeedback,
                context=context,
                is_child_algorithm=True,
            )
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            self.algRunner.runCreateSpatialIndex(
                otherRoadsLyr, context, multiStepFeedback, is_child_algorithm=True
            )
            currentStep += 1
        else:
            otherRoadsLyr = None

        multiStepFeedback.pushInfo(self.tr("Getting drainages"))
        multiStepFeedback.setCurrentStep(currentStep)
        drainagesWithNameLyr = self.algRunner.runFilterExpression(
            inputLyr=drainagesLyr,
            expression=f'"{drainageNameAttribute}" is not NULL',
            context=context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        currentStep += 1
        multiStepFeedback.setCurrentStep(currentStep)
        drainagesWithNameLyr = self.algRunner.runReprojectLayer(
            layer=drainagesWithNameLyr,
            targetCrs=QgsCoordinateReferenceSystem("EPSG:4674"),
            feedback=multiStepFeedback,
            context=context,
            is_child_algorithm=True,
        )

        currentStep += 1
        multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            drainagesWithNameLyr, context, multiStepFeedback, is_child_algorithm=True
        )
        currentStep += 1

        if getDrainagesWithoutName:
            multiStepFeedback.setCurrentStep(currentStep)
            drainagesWithoutNameLyr = self.algRunner.runFilterExpression(
                inputLyr=drainagesLyr,
                expression=f'"{drainageNameAttribute}" is NULL',
                context=context,
                feedback=multiStepFeedback,
                is_child_algorithm=True,
            )
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            drainagesWithoutNameLyr = self.algRunner.runReprojectLayer(
                layer=drainagesWithoutNameLyr,
                targetCrs=QgsCoordinateReferenceSystem("EPSG:4674"),
                feedback=multiStepFeedback,
                context=context,
                is_child_algorithm=True,
            )
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            self.algRunner.runCreateSpatialIndex(
                drainagesWithoutNameLyr,
                context,
                multiStepFeedback,
                is_child_algorithm=True,
            )
            currentStep += 1
        else:
            drainagesWithoutNameLyr = None

        for currentStepInsideLoop, feat in enumerate(
            geographicBoundaryLyr.getFeatures(), currentStep
        ):
            if multiStepFeedback.isCanceled():
                break
            multiStepFeedback.setCurrentStep(currentStepInsideLoop)
            self.currentStepText = self.tr(
                f"Evaluating region {currentStepInsideLoop-currentStep+1}/{nFeats}"
            )
            multiStepFeedback.setProgressText(self.currentStepText)
            localBoundsLyr = layerHandler.createMemoryLayerWithFeature(
                geographicBoundaryLyr, feat, context
            )
            featList = self.computePoints(
                contourLyr=parameters[self.CONTOUR_LINES],
                heightFieldName=heightFieldName,
                depressionExpression=depressionExpression,
                contourHeightInterval=contourHeightInterval,
                inputRaster=inputRaster,
                geographicBoundsLyr=localBoundsLyr,
                areaWithoutInformationLyr=areaWithoutInformationLyr
                if areaWithoutInformationLyr is not None
                else None,
                waterBodiesLyr=waterBodiesLyr,
                naturalPointFeaturesLyr=naturalPointFeaturesLyr
                if naturalPointFeaturesLyr is not None
                else None,
                drainagesWithNameLyr=drainagesWithNameLyr
                if drainagesWithNameLyr is not None
                else None,
                drainagesWithoutNameLyr=drainagesWithoutNameLyr
                if drainagesWithoutNameLyr is not None
                else None,
                mainRoadsLyr=mainRoadsLyr if mainRoadsLyr is not None else None,
                otherRoadsLyr=otherRoadsLyr if otherRoadsLyr is not None else None,
                onlyHilltops=onlyHilltops,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
            self.sink.addFeatures(featList, QgsFeatureSink.FastInsert)
        return {
            "OUTPUT": self.sink_id,
        }

    def computePoints(
        self,
        contourLyr,
        heightFieldName,
        contourHeightInterval,
        depressionExpression,
        inputRaster,
        geographicBoundsLyr,
        areaWithoutInformationLyr,
        waterBodiesLyr,
        naturalPointFeaturesLyr,
        drainagesWithNameLyr,
        drainagesWithoutNameLyr,
        mainRoadsLyr,
        otherRoadsLyr,
        onlyHilltops,
        fields,
        context,
        feedback,
    ):
        layerHandler = LayerHandler()
        nSteps = (
            17
            + (naturalPointFeaturesLyr is not None)
            + 3 * (waterBodiesLyr is not None)
            + (areaWithoutInformationLyr is not None)
        )  # handle this count after alg is done
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(nSteps, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(self.tr("Clipping raster"))
        frameCentroid = (
            [i for i in geographicBoundsLyr.getFeatures()][0].geometry().centroid()
        )
        originEpsg = QgsCoordinateReferenceSystem(
            geometryHandler.getSirgasAuthIdByPointLatLong(*frameCentroid.asPoint())
        )
        localGridSize = geometryHandler.convertDistance(
            self.gridSpacing,
            originEpsg=originEpsg,
            destinationEpsg=geographicBoundsLyr.crs(),
        )
        minusBufferedGeographicBoundsLyr = self.algRunner.runBuffer(
            inputLayer=geographicBoundsLyr,
            distance=-localGridSize / 64,
            endCapStyle=1,
            context=context,
        )
        clippedRaster = self.algRunner.runClipRasterLayer(
            inputRaster,
            mask=geographicBoundsLyr,
            context=context,
            feedback=feedback,
            nodata=-9999,
            outputRaster=QgsProcessingUtils.generateTempFilename(
                f"clip_{str(uuid4().hex)}.tif"
            ),
        )
        clippedRasterLyr = QgsProcessingUtils.mapLayerFromString(clippedRaster, context)
        localBufferDistance = geometryHandler.convertDistance(
            self.bufferDist,
            originEpsg=originEpsg,
            destinationEpsg=geographicBoundsLyr.crs(),
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(self.tr("Reading raster with numpy..."))
        localContourBufferLength = geometryHandler.convertDistance(
            self.contourBufferLength,
            originEpsg=originEpsg,
            destinationEpsg=geographicBoundsLyr.crs(),
        )
        contourBufferLyr = self.algRunner.runBuffer(
            inputLayer=contourLyr,
            distance=localContourBufferLength,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        if waterBodiesLyr is not None:
            localWaterBodiesLyr = self.algRunner.runExtractByLocation(
                inputLyr=waterBodiesLyr,
                intersectLyr=geographicBoundsLyr,
                context=context,
                feedback=multiStepFeedback,
                is_child_algorithm=True,
            )
            if multiStepFeedback is not None:
                currentStep += 1
                multiStepFeedback.setCurrentStep(currentStep)
            localBufferedWaterBodiesLyr = self.algRunner.runBuffer(
                inputLayer=localWaterBodiesLyr,
                distance=localContourBufferLength,
                context=context,
                feedback=multiStepFeedback,
                is_child_algorithm=True,
            )
            if multiStepFeedback is not None:
                currentStep += 1
                multiStepFeedback.setCurrentStep(currentStep)
            waterBodiesLyr = self.algRunner.runClip(
                localBufferedWaterBodiesLyr,
                overlayLayer=geographicBoundsLyr,
                context=context,
                feedback=multiStepFeedback,
            )
        if areaWithoutInformationLyr is not None:
            if multiStepFeedback is not None:
                currentStep += 1
                multiStepFeedback.setCurrentStep(currentStep)
                multiStepFeedback.pushInfo(
                    self.tr(
                        f"{self.currentStepText}: Running clip on area without information..."
                    )
                )
            areaWithoutInformationLyr = self.algRunner.runClip(
                areaWithoutInformationLyr,
                overlayLayer=geographicBoundsLyr,
                context=context,
                feedback=multiStepFeedback,
            )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr(f"{self.currentStepText}: Building masked raster...")
            )
        npRaster, maskedNpRaster, transform, maskLyr = self.readAndMaskRaster(
            clippedRasterLyr,
            geographicBoundsLyr,
            areaWithoutInformationLyr,
            waterBodiesLyr,
            context,
            self.algRunner,
            feedback=multiStepFeedback,
        )

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr("Getting max min feats and building exclusion polygons...")
            )
        maxFeats = self.getMaxFeatures(
            fields,
            npRaster,
            transform,
            distance=localBufferDistance,
            maskLyr=maskLyr,
            crs=maskLyr.crs() if maskLyr is not None else clippedRasterLyr.crs(),
            context=context,
            feedback=multiStepFeedback,
        )
        minFeats = self.getMinFeatures(
            fields,
            npRaster,
            transform,
            distance=localBufferDistance,
            maskLyr=maskLyr,
            crs=maskLyr.crs() if maskLyr is not None else clippedRasterLyr.crs(),
            context=context,
            feedback=multiStepFeedback,
        )
        maskLyr = (
            self.algRunner.runMergeVectorLayers([maskLyr, contourBufferLyr], context)
            if maskLyr is not None
            else contourBufferLyr
        )
        minMaxFeats = maxFeats + minFeats
        elevationPointsLayer = layerHandler.createMemoryLayerWithFeatures(
            featList=minMaxFeats,
            fields=fields,
            crs=clippedRasterLyr.crs(),
            wkbType=QgsWkbTypes.Point,
            context=context,
        )
        # compute number of points
        minNPoints, maxNPoints = self.getRangeOfNumberOfPoints(minMaxFeats)
        maxPointsPerGridUnit = maxNPoints // 9
        exclusionLyr = self.buildExclusionLyr(
            elevationPointsLayer,
            areaWithoutInformationLyr,
            waterBodiesLyr,
            localBufferDistance,
            context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(self.tr("Preparing contours..."))
        contourAreaDict, polygonLyr = self.prepareContours(
            contourLyr=contourLyr,
            geographicBoundsLyr=geographicBoundsLyr,
            heightFieldName=heightFieldName,
            depressionExpression=depressionExpression,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)

        if onlyHilltops:
            if multiStepFeedback is not None:
                multiStepFeedback.pushInfo(
                    self.tr("Extracting only hilltop points (all hilltops)...")
                )
            allHilltopPoints = self.extractAllHilltopPoints(
                contourLyr=contourLyr,
                polygonLyr=polygonLyr,
                rasterLyr=clippedRaster,
                geographicBoundsLyr=geographicBoundsLyr,
                minusBuffergeographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                originEpsg=originEpsg,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
            return allHilltopPoints

        # create grid
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        gridLyr, gridDict, gridLyrDict = self.buildGrid(
            geographicBoundsLyr, context, feedback=multiStepFeedback
        )

        # create points from first criteria
        if naturalPointFeaturesLyr is not None:
            if multiStepFeedback is not None:
                multiStepFeedback.pushInfo(
                    self.tr("Getting elevation points from natural points...")
                )
            elevationPointsFromNaturalPointFeatures = (
                self.getElevationPointsFromNaturalPoints(
                    npRaster=npRaster,
                    transform=transform,
                    naturalPointFeaturesLyr=naturalPointFeaturesLyr,
                    exclusionLyr=exclusionLyr,
                    polygonLyr=polygonLyr,
                    bufferDistance=localBufferDistance,
                    contourAreaDict=contourAreaDict,
                    contourHeightInterval=contourHeightInterval,
                    fields=fields,
                    gridLyr=gridLyr,
                    gridDict=gridDict,
                    gridLyrDict=gridLyrDict,
                    maxPointsPerGridUnit=maxPointsPerGridUnit,
                    context=context,
                    feedback=multiStepFeedback,
                )
            )

            self.updateExclusionLyr(
                exclusionLyr,
                elevationPointsFromNaturalPointFeatures,
                distance=localBufferDistance,
                context=context,
            )
            self.addPointsToMemoryLayer(
                elevationPointsLayer, elevationPointsFromNaturalPointFeatures, context
            )
            if multiStepFeedback is not None:
                currentStep += 1
                multiStepFeedback.setCurrentStep(currentStep)

        if multiStepFeedback is not None:
            multiStepFeedback.pushInfo(
                self.tr("Getting elevation points from hilltops...")
            )

        # create points from hilltops
        elevationPointsFromHilltops = self.getElevationPointsFromHilltops(
            rasterLyr=clippedRaster,
            contourLyr=contourLyr,
            polygonLyr=polygonLyr,
            exclusionLyr=exclusionLyr,
            geographicBoundsLyr=geographicBoundsLyr,
            depressionExpression=depressionExpression,
            minusBuffergeographicBoundsLyr=minusBufferedGeographicBoundsLyr,
            bufferDistance=localBufferDistance,
            originEpsg=originEpsg,
            contourAreaDict=contourAreaDict,
            contourHeightInterval=contourHeightInterval,
            gridLyr=gridLyr,
            gridDict=gridDict,
            maxPointsPerGridUnit=maxPointsPerGridUnit,
            fields=fields,
            context=context,
            feedback=multiStepFeedback,
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromHilltops,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer, elevationPointsFromHilltops, context
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr("Getting elevation points from main road intersections...")
            )
        elevationPointsFromRoadIntersections = (
            self.getElevationPointsFromLineIntersections(
                lineLyr1=mainRoadsLyr,
                lineLyr2=mainRoadsLyr,
                geographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                npRaster=npRaster,
                transform=transform,
                polygonLyr=polygonLyr,
                bufferDistance=localBufferDistance,
                exclusionLyr=exclusionLyr,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromRoadIntersections,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer, elevationPointsFromRoadIntersections, context
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr("Getting elevation points from other road intersections...")
            )
        elevationPointsFromOtherRoadIntersections = (
            self.getElevationPointsFromLineIntersections(
                lineLyr1=otherRoadsLyr,
                lineLyr2=otherRoadsLyr,
                geographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                npRaster=npRaster,
                transform=transform,
                polygonLyr=polygonLyr,
                bufferDistance=localBufferDistance,
                exclusionLyr=exclusionLyr,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromOtherRoadIntersections,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer, elevationPointsFromOtherRoadIntersections, context
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr(
                    "Getting elevation points from intersections of road and rivers with names outside polygons..."
                )
            )
        elevationPointsFromIntersectionsOfRiverWithNamesAndRoads = (
            self.getElevationPointsFromLineIntersections(
                lineLyr1=drainagesWithNameLyr,
                lineLyr2=mainRoadsLyr,
                geographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                npRaster=npRaster,
                transform=transform,
                polygonLyr=polygonLyr,
                bufferDistance=localBufferDistance,
                exclusionLyr=exclusionLyr,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromIntersectionsOfRiverWithNamesAndRoads,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer,
            elevationPointsFromIntersectionsOfRiverWithNamesAndRoads,
            context,
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr(
                    "Getting elevation points from intersections of road and rivers without names..."
                )
            )

        elevationPointsFromIntersectionsOfRiversWithoutNamesAndRoads = (
            self.getElevationPointsFromLineIntersections(
                lineLyr1=drainagesWithoutNameLyr,
                lineLyr2=mainRoadsLyr,
                geographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                npRaster=npRaster,
                transform=transform,
                polygonLyr=polygonLyr,
                bufferDistance=localBufferDistance,
                exclusionLyr=exclusionLyr,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromIntersectionsOfRiversWithoutNamesAndRoads,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer,
            elevationPointsFromIntersectionsOfRiversWithoutNamesAndRoads,
            context,
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr(
                    "Getting elevation points from intersections of rivers with names..."
                )
            )

        elevationPointsFromIntersectionsOfRiversWithNames = (
            self.getElevationPointsFromLineIntersections(
                lineLyr1=drainagesWithNameLyr,
                lineLyr2=drainagesWithNameLyr,
                geographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                npRaster=npRaster,
                transform=transform,
                polygonLyr=polygonLyr,
                bufferDistance=localBufferDistance,
                exclusionLyr=exclusionLyr,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromIntersectionsOfRiversWithNames,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer,
            elevationPointsFromIntersectionsOfRiversWithNames,
            context,
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr(
                    "Getting elevation points from intersections of rivers with names..."
                )
            )

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr(
                    "Getting elevation points from intersections of rivers with and without names..."
                )
            )

        elevationPointsFromIntersectionsOfRiversWithAndWithoutNames = (
            self.getElevationPointsFromLineIntersections(
                lineLyr1=drainagesWithNameLyr,
                lineLyr2=drainagesWithoutNameLyr,
                geographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                npRaster=npRaster,
                transform=transform,
                polygonLyr=polygonLyr,
                bufferDistance=localBufferDistance,
                exclusionLyr=exclusionLyr,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromIntersectionsOfRiversWithAndWithoutNames,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer,
            elevationPointsFromIntersectionsOfRiversWithAndWithoutNames,
            context,
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr(
                    "Getting elevation points from intersections of rivers without names..."
                )
            )

        elevationPointsFromIntersectionsOfRiversWithoutNames = (
            self.getElevationPointsFromLineIntersections(
                lineLyr1=drainagesWithoutNameLyr,
                lineLyr2=drainagesWithoutNameLyr,
                geographicBoundsLyr=minusBufferedGeographicBoundsLyr,
                npRaster=npRaster,
                transform=transform,
                polygonLyr=polygonLyr,
                bufferDistance=localBufferDistance,
                exclusionLyr=exclusionLyr,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
        )
        self.updateExclusionLyr(
            exclusionLyr,
            elevationPointsFromIntersectionsOfRiversWithoutNames,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer,
            elevationPointsFromIntersectionsOfRiversWithoutNames,
            context,
        )
        if elevationPointsLayer.featureCount() > maxNPoints:
            return elevationPointsLayer.getFeatures()

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(
                self.tr("Getting elevation points from plane areas...")
            )

        planeAreasElevationPoints = self.getElevationPointsFromPlaneAreas(
            contourLyr=contourLyr,
            geographicBoundsLyr=geographicBoundsLyr,
            rasterLyr=clippedRasterLyr,
            exclusionLyr=exclusionLyr,
            polygonLyr=polygonLyr,
            originEpsg=originEpsg,
            bufferDistance=localBufferDistance,
            contourAreaDict=contourAreaDict,
            contourHeightInterval=contourHeightInterval,
            gridLyr=gridLyr,
            gridDict=gridDict,
            maxPointsPerGridUnit=maxPointsPerGridUnit,
            fields=fields,
            context=context,
            feedback=multiStepFeedback,
        )
        self.updateExclusionLyr(
            exclusionLyr,
            planeAreasElevationPoints,
            distance=localBufferDistance,
            context=context,
        )
        self.addPointsToMemoryLayer(
            elevationPointsLayer, planeAreasElevationPoints, context
        )
        return elevationPointsLayer.getFeatures()

    def readAndMaskRaster(
        self,
        clippedRasterLyr,
        geographicBoundsLyr,
        areaWithoutInformationLyr,
        waterBodiesLyr,
        context,
        algRunner,
        feedback,
    ):
        nSteps = 6
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(nSteps, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        ds, npRaster = rasterHandler.readAsNumpy(clippedRasterLyr)
        transform = rasterHandler.getCoordinateTransform(ds)
        nanIndexes = np.isnan(npRaster)
        npRaster = (np.rint(npRaster)).astype(float)
        npRaster[nanIndexes] = np.nan
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        layerList = list(
            filter(
                lambda x: x is not None,
                [areaWithoutInformationLyr, waterBodiesLyr],
            )
        )
        if layerList == []:
            return npRaster, npRaster, transform, None
        maskLyr = (
            self.algRunner.runMergeVectorLayers(
                layerList, context, feedback=multiStepFeedback
            )
            if len(layerList) > 1
            else layerList[0]
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        maskLyr = self.algRunner.runClip(
            maskLyr, geographicBoundsLyr, context, feedback=multiStepFeedback
        )

        if maskLyr.featureCount() == 0:
            return npRaster, npRaster, transform, None
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        maskedNpRaster = self.maskFeaturesFromLayerOnRaster(
            rasterLyr=clippedRasterLyr,
            geographicBoundsLyr=geographicBoundsLyr,
            maskLyr=maskLyr,
            context=context,
            feedback=multiStepFeedback,
            npRaster=npRaster,
        )
        nanIndexes = np.isnan(maskedNpRaster)
        maskedNpRaster = (np.rint(maskedNpRaster)).astype(float)
        maskedNpRaster[nanIndexes] = np.nan
        return npRaster, maskedNpRaster, transform, maskLyr

    def getRangeOfNumberOfPoints(self, featList):
        minValue = min(featList, key=lambda x: x["cota"])["cota"]
        maxValue = max(featList, key=lambda x: x["cota"])["cota"]
        if abs(maxValue - minValue) > 250:
            return (50, 150)
        return (25, 75)

    def buildExclusionLyr(
        self,
        elevationPointsLayer: QgsVectorLayer,
        areaWithoutInformationLyr: QgsVectorLayer,
        waterBodiesLyr: QgsVectorLayer,
        bufferDistance: float,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> QgsVectorLayer:
        buffer = self.algRunner.runBuffer(
            elevationPointsLayer,
            distance=bufferDistance,
            endCapStyle=0,
            joinStyle=0,
            context=context,
            feedback=feedback,
        )
        layerList = [buffer]
        for lyr in [areaWithoutInformationLyr, waterBodiesLyr]:
            if lyr is not None:
                auxLyr = self.algRunner.runMultipartToSingleParts(
                    lyr, context, is_child_algorithm=True
                )
                layerList.append(auxLyr)
        outputLyr = (
            buffer
            if len(layerList) == 1
            else self.algRunner.runMergeVectorLayers(layerList, context)
        )
        self.algRunner.runCreateSpatialIndex(
            outputLyr, context, is_child_algorithm=True
        )
        return outputLyr

    def addPointsToMemoryLayer(
        self, lyr: QgsVectorLayer, featList: List, context: QgsProcessingContext
    ) -> None:
        lyr.startEditing()
        lyr.beginEditCommand("adding features")
        lyr.addFeatures(featList)
        lyr.endEditCommand()
        lyr.commitChanges()
        self.algRunner.runCreateSpatialIndex(
            inputLyr=lyr, context=context, is_child_algorithm=True
        )

    def prepareContours(
        self,
        contourLyr: QgsVectorLayer,
        geographicBoundsLyr: QgsVectorLayer,
        heightFieldName: str,
        depressionExpression: str,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> Tuple[Dict, QgsVectorLayer]:
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(5, feedback)
            if feedback is not None
            else None
        )
        spatialRelatonsHandler = SpatialRelationsHandler()
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        clippedContourLyr = self.algRunner.runClip(
            inputLayer=contourLyr,
            overlayLayer=geographicBoundsLyr,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        clippedContourLyr = self.algRunner.runMultipartToSingleParts(
            inputLayer=clippedContourLyr,
            context=context,
            feedback=multiStepFeedback,
        )
        emptyList = [
            feat.id()
            for feat in clippedContourLyr.getFeatures()
            if feat.geometry().isEmpty()
        ]
        if emptyList != []:
            clippedContourLyr.startEditing()
            clippedContourLyr.beginEditCommand("deleting empty")
            clippedContourLyr.deleteFeatures(emptyList)
            clippedContourLyr.endEditCommand()
            clippedContourLyr.commitChanges()
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        (
            contourSpatialIdx,
            contourIdDict,
            _,
            __,
        ) = spatialRelatonsHandler.buildSpatialIndexAndIdDictRelateNodesAndAttributeGroupDict(
            inputLyr=clippedContourLyr,
            attributeName=heightFieldName,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        contourAreaDict, polygonLyr = spatialRelatonsHandler.buildContourAreaDict(
            inputLyr=clippedContourLyr,
            geoBoundsLyr=geographicBoundsLyr,
            attributeName=heightFieldName,
            contourSpatialIdx=contourSpatialIdx,
            contourIdDict=contourIdDict,
            depressionExpression=depressionExpression,  # TODO
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            polygonLyr, context, feedback=multiStepFeedback, is_child_algorithm=True
        )
        return contourAreaDict, polygonLyr

    def buildGrid(
        self,
        geographicBoundsLyr: QgsVectorLayer,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> Tuple[QgsVectorLayer, Dict[QByteArray, int], Dict[QByteArray, QgsVectorLayer]]:
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(4, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        boundsGeom = [i for i in geographicBoundsLyr.getFeatures()][0].geometry()
        xmin, ymin, xmax, ymax = boundsGeom.boundingBox().toRectF().getCoords()
        hSpacing = abs(xmax - xmin) / 3
        vSpacing = abs(ymax - ymin) / 3
        gridLyr = self.algRunner.runCreateGrid(
            extent=geographicBoundsLyr.extent(),
            crs=geographicBoundsLyr.crs(),
            hSpacing=hSpacing,
            vSpacing=vSpacing,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        stepSize = 100 / 9
        gridDict = defaultdict(int)
        gridLyrDict = dict()
        for current, feat in enumerate(gridLyr.getFeatures()):
            if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                break
            geom = feat.geometry()
            geomWkb = geom.asWkb()
            gridDict[geomWkb] = 0
            gridLyrDict[geomWkb] = LayerHandler().createMemoryLayerWithFeature(
                geographicBoundsLyr, feat, context
            )
            if multiStepFeedback is not None:
                multiStepFeedback.setProgress(current * stepSize)
        return gridLyr, gridDict, gridLyrDict

    def getElevationPointsFromNaturalPoints(
        self,
        npRaster: np.array,
        transform: Affine,
        naturalPointFeaturesLyr: QgsVectorLayer,
        exclusionLyr: QgsVectorLayer,
        polygonLyr: QgsVectorLayer,
        bufferDistance: float,
        contourAreaDict: Dict,
        contourHeightInterval: float,
        gridLyr: QgsVectorLayer,
        gridDict: Dict[QByteArray, int],
        maxPointsPerGridUnit: int,
        fields: QgsFields,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> List[QgsFeature]:
        nFeats = naturalPointFeaturesLyr.featureCount()
        if nFeats == 0:
            return
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(2, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        pointList = self.getElevationPointsFromLayer(
            npRaster, transform, naturalPointFeaturesLyr, fields
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        return self.filterWithAllCriteria(
            inputPointList=pointList,
            referenceLyr=naturalPointFeaturesLyr,
            exclusionLyr=exclusionLyr,
            polygonLyr=polygonLyr,
            bufferDistance=bufferDistance,
            contourAreaDict=contourAreaDict,
            contourHeightInterval=contourHeightInterval,
            gridLyr=gridLyr,
            gridDict=gridDict,
            maxPointsPerGridUnit=maxPointsPerGridUnit,
            fields=fields,
            context=context,
            feedback=multiStepFeedback,
        )

    def getElevationPointsFromLayer(self, npRaster, transform, lyr, fields):
        pointList = rasterHandler.createFeatureListWithPointList(
            pointList=lyr.getFeatures(),
            fieldName="cota",
            fields=fields,
            npRaster=npRaster,
            transform=transform,
            defaultAtributeMap=dict(self.defaultAttrMap),
        )
        return filter(lambda x: x is not None, pointList)

    def validateGeneratedPointsAgainstTerrainModel(
        self,
        elevationPoints: List[QgsFeature],
        polygonLyr: QgsVectorLayer,
        contourAreaDict: Dict,
        threshold: float,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ):
        flagDict = SpatialRelationsHandler().findElevationPointsOutOfThreshold(
            elevationPoints=elevationPoints,
            polygonLyr=polygonLyr,
            contourAreaDict=contourAreaDict,
            threshold=threshold,
            elevationPointHeightFieldName="cota",
            context=context,
            feedback=feedback,
        )
        return (
            elevationPoints
            if flagDict == dict()
            else [
                point
                for point in elevationPoints
                if point.geometry().asWkb() not in flagDict
            ]
        )

    def getElevationPointsFromLineIntersections(
        self,
        lineLyr1,
        lineLyr2,
        geographicBoundsLyr,
        npRaster,
        transform,
        exclusionLyr: QgsVectorLayer,
        polygonLyr: QgsVectorLayer,
        bufferDistance: float,
        contourAreaDict: Dict,
        contourHeightInterval: float,
        gridLyr: QgsVectorLayer,
        gridDict: Dict[QByteArray, int],
        maxPointsPerGridUnit: int,
        fields: QgsFields,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ):
        if lineLyr1 is None or lineLyr2 is None:
            return []
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(17, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        localLineLyr1 = self.algRunner.runClip(
            lineLyr1,
            overlayLayer=geographicBoundsLyr,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            localLineLyr1, context, feedback=multiStepFeedback, is_child_algorithm=True
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        localLineLyr2 = self.algRunner.runClip(
            lineLyr2,
            overlayLayer=geographicBoundsLyr,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            localLineLyr2, context, feedback=multiStepFeedback, is_child_algorithm=True
        )
        if localLineLyr1.featureCount() == 0 or localLineLyr2.featureCount() == 0:
            return []
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        lineIntersectionPointsLyr = self.algRunner.runLineIntersections(
            localLineLyr1,
            intersectLyr=localLineLyr2,
            context=context,
            feedback=multiStepFeedback,
        )
        if lineIntersectionPointsLyr.featureCount() == 0:
            return []
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        intersectionBuffers = self.algRunner.runBuffer(
            inputLayer=lineIntersectionPointsLyr,
            distance=1e-5,
            context=context,
            is_child_algorithm=True,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            intersectionBuffers,
            context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        extractedLines1 = self.algRunner.runExtractByLocation(
            inputLyr=localLineLyr1,
            intersectLyr=intersectionBuffers,
            context=context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            extractedLines1,
            context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        extractedLines2 = self.algRunner.runExtractByLocation(
            inputLyr=localLineLyr2,
            intersectLyr=intersectionBuffers,
            context=context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            extractedLines2,
            context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        merged = self.algRunner.runMergeVectorLayers(
            inputList=[extractedLines1, extractedLines2],
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        intersectedLines = self.algRunner.runSplitLinesWithLines(
            inputLyr=merged,
            linesLyr=merged,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runRemoveDuplicatedGeometries(
            inputLyr=intersectedLines,
            context=context,
            feedback=multiStepFeedback,
        )
        intersectedLines.commitChanges()
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            intersectedLines,
            context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        filteredPointsIntersections = self.keepThirdOrderOrHigherPoints(
            points=lineIntersectionPointsLyr,
            linesLyr=intersectedLines,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        localAttrDefaultMap = dict(self.defaultAttrMap)
        localAttrDefaultMap["suprimir_simbologia"] = 1
        pointList = rasterHandler.createFeatureListWithPointList(
            pointList=filteredPointsIntersections,
            fieldName="cota",
            fields=fields,
            npRaster=npRaster,
            transform=transform,
            defaultAtributeMap=localAttrDefaultMap,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        return self.filterWithAllCriteria(
            inputPointList=filter(lambda x: x is not None, pointList),
            referenceLyr=localLineLyr1,
            exclusionLyr=exclusionLyr,
            polygonLyr=polygonLyr,
            bufferDistance=bufferDistance,
            contourAreaDict=contourAreaDict,
            contourHeightInterval=contourHeightInterval,
            gridLyr=gridLyr,
            gridDict=gridDict,
            maxPointsPerGridUnit=maxPointsPerGridUnit,
            fields=fields,
            context=context,
            feedback=multiStepFeedback,
        )

    def getMaxFeatures(
        self,
        fields,
        npRaster,
        transform,
        distance,
        maskLyr,
        crs,
        context,
        feedback=None,
    ):
        featSet = set()
        maxCoordinatesArray = rasterHandler.getMaxCoordinatesFromNpArray(npRaster)
        npRasterCopy = np.array(npRaster)
        if feedback is not None:
            feedback.pushInfo(
                self.tr("Creating max feature list from pixel coordinates array...")
            )
        maxFeatList = (
            rasterHandler.createFeatureListWithPixelValuesFromPixelCoordinatesArray(
                maxCoordinatesArray,
                fieldName="cota",
                fields=fields,
                npRaster=npRaster,
                transform=transform,
                defaultAtributeMap=dict(self.defaultAttrMap),
            )
        )
        cotaMax = min(f["cota"] for f in maxFeatList)
        while True:
            maxFeatLyr = LayerHandler().createMemoryLayerWithFeatures(
                featList=maxFeatList,
                fields=fields,
                crs=crs,
                wkbType=QgsWkbTypes.Point,
                context=QgsProcessingContext(),
            )
            filteredFeatureList = self.filterFeaturesByDistanceAndExclusionLayer(
                candidatesPointLyr=maxFeatLyr,
                exclusionLyr=maskLyr,
                distance=distance / 10,
                context=QgsProcessingContext(),
                feedback=None,
            )
            if filteredFeatureList != []:
                featSet |= self.filterFeaturesByBuffer(
                    filteredFeatureList, distance / 10, cotaMaisAlta=True
                )
                break
            if len(npRasterCopy[npRasterCopy == cotaMax]) > 0:
                npRasterCopy[npRasterCopy == cotaMax] = np.nan
            maxCoordinatesArray = rasterHandler.getMaxCoordinatesFromNpArray(
                npRasterCopy
            )
            if maxCoordinatesArray.size == 0:
                return list(featSet)
            maxFeatList = (
                rasterHandler.createFeatureListWithPixelValuesFromPixelCoordinatesArray(
                    maxCoordinatesArray,
                    fieldName="cota",
                    fields=fields,
                    npRaster=npRaster,
                    transform=transform,
                    defaultAtributeMap=dict(self.defaultAttrMap),
                )
            )
            cotaMax -= 1
            if feedback is not None and feedback.isCanceled():
                break
        candidatesPointLyr = LayerHandler().createMemoryLayerWithFeatures(
            featList=list(featSet),
            fields=fields,
            crs=crs,
            wkbType=QgsWkbTypes.Point,
            context=context,
        )
        filteredPoints = self.filterFeaturesByDistanceAndExclusionLayer(
            candidatesPointLyr=candidatesPointLyr,
            exclusionLyr=maskLyr,
            distance=10 * distance,
            context=context,
            feedback=None,
        )
        return filteredPoints

    def getMinFeatures(
        self,
        fields,
        npRaster,
        transform,
        distance,
        maskLyr,
        crs,
        context,
        feedback=None,
    ):
        featSet = set()
        npRasterCopy = np.array(npRaster)
        minCoordinatesArray = rasterHandler.getMinCoordinatesFromNpArray(npRasterCopy)
        if len(minCoordinatesArray) == 0:
            return list(featSet)
        if feedback is not None:
            feedback.pushInfo(
                self.tr("Creating min feature list from pixel coordinates array...")
            )
        minFeatList = (
            rasterHandler.createFeatureListWithPixelValuesFromPixelCoordinatesArray(
                minCoordinatesArray,
                fieldName="cota",
                fields=fields,
                npRaster=npRasterCopy,
                transform=transform,
                defaultAtributeMap=dict(self.defaultAttrMap),
            )
        )
        cota = min(f["cota"] for f in minFeatList)
        while True:
            minFeatLyr = LayerHandler().createMemoryLayerWithFeatures(
                featList=minFeatList,
                fields=fields,
                crs=crs,
                wkbType=QgsWkbTypes.Point,
                context=QgsProcessingContext(),
            )
            filteredFeatureList = self.filterFeaturesByDistanceAndExclusionLayer(
                candidatesPointLyr=minFeatLyr,
                exclusionLyr=maskLyr,
                distance=distance,
                context=QgsProcessingContext(),
                feedback=None,
            )
            if filteredFeatureList != []:
                featSet |= set(filteredFeatureList)
                break
            npRasterCopy[npRasterCopy == cota] = np.nan
            minCoordinatesArray = rasterHandler.getMinCoordinatesFromNpArray(
                npRasterCopy
            )
            if minCoordinatesArray.size == 0:
                break
            minFeatList = (
                rasterHandler.createFeatureListWithPixelValuesFromPixelCoordinatesArray(
                    minCoordinatesArray,
                    fieldName="cota",
                    fields=fields,
                    npRaster=npRasterCopy,
                    transform=transform,
                    defaultAtributeMap=dict(self.defaultAttrMap),
                )
            )
            cota += 1
            if feedback is not None and feedback.isCanceled():
                break
        candidatesPointLyr = LayerHandler().createMemoryLayerWithFeatures(
            featList=minFeatList,
            fields=fields,
            crs=crs,
            wkbType=QgsWkbTypes.Point,
            context=context,
        )
        filteredPoints = self.filterFeaturesByDistanceAndExclusionLayer(
            candidatesPointLyr=candidatesPointLyr,
            exclusionLyr=maskLyr,
            distance=10 * distance,
            context=context,
            feedback=None,
        )
        return filteredPoints

    def filterWithAllCriteria(
        self,
        inputPointList: List[QgsFeature],
        referenceLyr: QgsVectorLayer,
        exclusionLyr: QgsVectorLayer,
        polygonLyr: QgsVectorLayer,
        bufferDistance: float,
        contourAreaDict: Dict,
        contourHeightInterval: float,
        gridLyr: QgsVectorLayer,
        gridDict: Dict[QByteArray, int],
        maxPointsPerGridUnit: int,
        fields: QgsFields,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ):
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(3, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        candidatesPointLyr = LayerHandler().createMemoryLayerWithFeatures(
            featList=inputPointList,
            fields=fields,
            crs=referenceLyr.crs(),
            wkbType=QgsWkbTypes.Point,
            context=context,
        )
        pointList = self.filterFeaturesByDistanceAndExclusionLayer(
            candidatesPointLyr=candidatesPointLyr,
            exclusionLyr=exclusionLyr,
            distance=bufferDistance,
            context=context,
            feedback=multiStepFeedback,
        )
        if len(pointList) == []:
            return []
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        pointList = self.validateGeneratedPointsAgainstTerrainModel(
            elevationPoints=pointList,
            polygonLyr=polygonLyr,
            contourAreaDict=contourAreaDict,
            threshold=contourHeightInterval,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        newCandidatesPointLyr = LayerHandler().createMemoryLayerWithFeatures(
            featList=pointList,
            fields=fields,
            crs=referenceLyr.crs(),
            wkbType=QgsWkbTypes.Point,
            context=context,
        )
        return self.filterPointsAgainstGrid(
            pointLyr=newCandidatesPointLyr,
            gridLyr=gridLyr,
            gridDict=gridDict,
            maxPointsPerGridUnit=maxPointsPerGridUnit,
            feedback=multiStepFeedback,
        )

    def filterFeaturesByBuffer(
        self, filterFeatList: List, distance, cotaMaisAlta=False
    ):
        outputSet = set()
        exclusionGeom = None
        for feat in filterFeatList:
            geom = feat.geometry()
            buffer = geom.buffer(
                distance, 10, Qgis.EndCapStyle.Round, Qgis.JoinStyle.Round, -1
            )
            if exclusionGeom is not None and exclusionGeom.intersects(geom):
                continue
            feat["cota_mais_alta"] = 1 if cotaMaisAlta else 2
            exclusionGeom = (
                buffer if exclusionGeom is None else exclusionGeom.combine(buffer)
            )
            outputSet.add(feat)
        return outputSet

    def filterFeaturesByDistanceAndExclusionLayer(
        self,
        candidatesPointLyr: QgsVectorLayer,
        exclusionLyr: QgsVectorLayer,
        distance: float,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> List:
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(5, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)

        candidatesPointLyrWithId = self.algRunner.runCreateFieldWithExpression(
            inputLyr=candidatesPointLyr,
            expression="$id",
            fieldName="featid",
            fieldType=1,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            candidatesPointLyrWithId,
            context,
            feedback=multiStepFeedback,
            is_child_algorithm=True,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        disjointLyr = (
            self.algRunner.runExtractByLocation(
                candidatesPointLyrWithId,
                exclusionLyr,
                predicate=[AlgRunner.Disjoint],
                context=context,
                feedback=multiStepFeedback,
            )
            if exclusionLyr is not None
            else candidatesPointLyrWithId
        )
        nFeats = disjointLyr.featureCount()
        if nFeats == 0:
            return []
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        self.algRunner.runCreateSpatialIndex(
            disjointLyr, context, feedback=multiStepFeedback, is_child_algorithm=True
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        outputSet, exclusionSet = set(), set()
        stepSize = 100 / nFeats
        request = QgsFeatureRequest()
        orderByClause1 = QgsFeatureRequest.OrderByClause(
            "cota_mais_alta", ascending=True
        )
        orderByClause2 = QgsFeatureRequest.OrderByClause("cota", ascending=False)
        # ordenamento geográfico, da esquerda para direita, de cima para baixo
        orderByClause3 = QgsFeatureRequest.OrderByClause("$y", ascending=False)
        orderByClause4 = QgsFeatureRequest.OrderByClause("$x", ascending=True)
        orderby = QgsFeatureRequest.OrderBy(
            [orderByClause1, orderByClause2, orderByClause3, orderByClause4]
        )
        request.setOrderBy(orderby)
        for current, feat in enumerate(disjointLyr.getFeatures(request)):
            if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                break
            if feat["featid"] in exclusionSet:
                if multiStepFeedback is not None:
                    multiStepFeedback.setProgress(current * stepSize)
                continue
            outputSet.add(feat)
            geom = feat.geometry()
            buffer = geom.buffer(
                distance, 10, Qgis.EndCapStyle.Round, Qgis.JoinStyle.Round, -1
            )
            bbox = buffer.boundingBox()
            for candidateFeat in disjointLyr.getFeatures(bbox):
                if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                    break
                if candidateFeat["featid"] == feat["featid"]:
                    continue
                if not candidateFeat.geometry().intersects(buffer):
                    continue
                exclusionSet.add(candidateFeat["featid"])
            if multiStepFeedback is not None:
                multiStepFeedback.setProgress(current * stepSize)
        return list(outputSet)

    def filterPointsAgainstGrid(
        self,
        pointLyr: QgsVectorLayer,
        gridLyr: QgsVectorLayer,
        gridDict: Dict[QByteArray, int],
        maxPointsPerGridUnit: int,
        feedback: QgsFeedback,
    ) -> List[QgsFeature]:
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(9, feedback)
            if feedback is not None
            else None
        )
        outputSet = set()
        for currentStep, feat in enumerate(gridLyr.getFeatures()):
            if multiStepFeedback is not None:
                multiStepFeedback.setCurrentStep(currentStep)
            if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                break
            geom = feat.geometry()
            bbox = geom.boundingBox()
            geomWkb = geom.asWkb()
            if gridDict[geomWkb] >= maxPointsPerGridUnit:
                continue
            gridFeats = [
                f for f in pointLyr.getFeatures(bbox) if f.geometry().intersects(geom)
            ]
            if gridFeats == []:
                continue
            selectedGridFeatsSet = set(
                islice(gridFeats, maxPointsPerGridUnit - gridDict[geomWkb])
            )
            gridDict[geomWkb] += len(selectedGridFeatsSet)
            outputSet |= selectedGridFeatsSet
        return list(outputSet)

    def updateExclusionLyr(
        self,
        exclusionLyr: QgsVectorLayer,
        pointList: List,
        distance: float,
        context: QgsProcessingContext,
    ) -> None:
        exclusionLyr.startEditing()
        exclusionLyr.beginEditCommand("updating exclusion layer")
        for feat in pointList:
            geom = feat.geometry()
            buffer = (
                geom.buffer(
                    distance, 10, Qgis.EndCapStyle.Round, Qgis.JoinStyle.Round, -1
                )
                if geom.type() == QgsWkbTypes.GeometryType.PointGeometry
                else geom
            )
            newFeat = QgsVectorLayerUtils.createFeature(exclusionLyr, buffer)
            exclusionLyr.addFeature(newFeat)
        exclusionLyr.endEditCommand()
        exclusionLyr.commitChanges()
        self.algRunner.runCreateSpatialIndex(
            exclusionLyr, context, is_child_algorithm=True
        )

    def maskFeaturesFromLayerOnRaster(
        self,
        rasterLyr,
        geographicBoundsLyr,
        maskLyr,
        context,
        feedback,
        npRaster,
    ):
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(2, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            multiStepFeedback.setCurrentStep(0)
        clippedArea = self.algRunner.runClip(
            maskLyr,
            geographicBoundsLyr,
            context,
            multiStepFeedback,
        )
        if multiStepFeedback is not None:
            multiStepFeedback.setCurrentStep(1)
        npMask = (
            None
            if clippedArea.featureCount() == 0
            else rasterHandler.buildNumpyNodataMask(
                rasterLyr=rasterLyr, vectorLyr=clippedArea
            )
        )
        if npMask is not None:
            npMask.resize(npRaster.shape, refcheck=False)
            npRaster = npRaster + npMask
        return npRaster

    def getElevationPointsFromHilltops(
        self,
        rasterLyr: QgsRasterLayer,
        contourLyr: QgsVectorLayer,
        polygonLyr: QgsVectorLayer,
        exclusionLyr: QgsVectorLayer,
        geographicBoundsLyr: QgsVectorLayer,
        depressionExpression: str,
        minusBuffergeographicBoundsLyr: QgsVectorLayer,
        bufferDistance: float,
        originEpsg: QgsCoordinateReferenceSystem,
        contourAreaDict: Dict,
        contourHeightInterval: float,
        gridLyr: QgsVectorLayer,
        gridDict: Dict[QByteArray, int],
        maxPointsPerGridUnit: int,
        fields: QgsFields,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> List[QgsFeature]:
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(2, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        pointList, hillTopsLyr = self.extractElevationPointsFromHilltops(
            contourLayer=contourLyr,
            polygonLayer=polygonLyr,
            rasterLyr=rasterLyr,
            geographicBoundsLyr=geographicBoundsLyr,
            minusBuffergeographicBoundsLyr=minusBuffergeographicBoundsLyr,
            originEpsg=originEpsg,
            fields=fields,
            context=context,
            feedback=multiStepFeedback,
        )

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        depressionHilltops = [
            feat for feat in hillTopsLyr.getFeatures(depressionExpression)
        ]
        if len(depressionHilltops) > 0:
            depressionPoints = [
                f
                for f in pointList
                if any(
                    f.geometry().intersects(i.geometry()) for i in depressionHilltops
                )
            ]
            filteredDepressionPoints = self.filterWithAllCriteria(
                inputPointList=depressionPoints,
                referenceLyr=polygonLyr,
                exclusionLyr=exclusionLyr,
                polygonLyr=polygonLyr,
                bufferDistance=bufferDistance,
                contourAreaDict=contourAreaDict,
                contourHeightInterval=contourHeightInterval,
                gridLyr=gridLyr,
                gridDict=gridDict,
                maxPointsPerGridUnit=maxPointsPerGridUnit,
                fields=fields,
                context=context,
                feedback=multiStepFeedback,
            )
            self.updateExclusionLyr(
                exclusionLyr=exclusionLyr,
                pointList=filteredDepressionPoints,
                distance=bufferDistance,
                context=context,
            )
        else:
            filteredDepressionPoints = []
        filteredPoints = self.filterWithAllCriteria(
            inputPointList=list(set(pointList) - set(filteredDepressionPoints)),
            referenceLyr=polygonLyr,
            exclusionLyr=exclusionLyr,
            polygonLyr=polygonLyr,
            bufferDistance=bufferDistance,
            contourAreaDict=contourAreaDict,
            contourHeightInterval=contourHeightInterval,
            gridLyr=gridLyr,
            gridDict=gridDict,
            maxPointsPerGridUnit=maxPointsPerGridUnit,
            fields=fields,
            context=context,
            feedback=multiStepFeedback,
        )
        return filteredDepressionPoints + filteredPoints

    def extractElevationPointsFromHilltops(
        self,
        contourLayer: QgsVectorLayer,
        polygonLayer: QgsVectorLayer,
        geographicBoundsLyr: QgsVectorLayer,
        minusBuffergeographicBoundsLyr: QgsVectorLayer,
        rasterLyr: QgsRasterLayer,
        originEpsg: QgsCoordinateReferenceSystem,
        fields: QgsFields,
        context: QgsProcessingContext,
        feedback: QgsFeedback = None,
    ) -> List:
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(3, feedback)
            if feedback is not None
            else None
        )
        spatialRelationsHandler = SpatialRelationsHandler()
        layerHandler = LayerHandler()
        currentStep = 0
        if multiStepFeedback is not None:
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(self.tr("Finding hilltops..."))
        hillTopsLyr = spatialRelationsHandler.createHilltopLayerFromPolygonLayer(
            polygonLayer=polygonLayer,
            geographicBoundsLyr=geographicBoundsLyr,
            context=context,
            feedback=multiStepFeedback,
            computeOrder=True,
        )
        self.algRunner.runCreateSpatialIndex(
            hillTopsLyr, context=context, is_child_algorithm=True
        )
        if multiStepFeedback is not None:
            multiStepFeedback.setCurrentStep(currentStep)
            currentStep += 1
        minContourLength = geometryHandler.convertDistance(
            self.minContourLength,
            originEpsg=originEpsg,
            destinationEpsg=geographicBoundsLyr.crs(),
        )
        minusBufferLength = geometryHandler.convertDistance(
            self.contourBufferLength,
            originEpsg=originEpsg,
            destinationEpsg=geographicBoundsLyr.crs(),
        )
        clippedHilltops = self.algRunner.runClip(
            inputLayer=hillTopsLyr,
            overlayLayer=minusBuffergeographicBoundsLyr,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        nFeats = clippedHilltops.featureCount()
        if nFeats == 0:
            return [], clippedHilltops
        self.algRunner.runCreateSpatialIndex(
            clippedHilltops, context=context, is_child_algorithm=True
        )
        clippedHilltops = self.algRunner.runJoinAttributesByLocation(
            inputLyr=clippedHilltops,
            joinLyr=contourLayer,
            context=context,
            discardNonMatching=True,
            method=1,
        )
        stepSize = 100 / nFeats
        featList = []
        minArea = minContourLength**2 / (4 * math.pi)
        hilltopList = sorted(
            filter(
                lambda x: x["order_count"] > 1
                and x.geometry().length() > minContourLength
                and x.geometry().area() > minArea,
                clippedHilltops.getFeatures(),
            ),
            key=lambda x: (x["order_count"], 1.0 / x.geometry().area()),
            reverse=True,
        )
        for current, hilltopFeat in enumerate(hilltopList):
            if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                break
            geom = hilltopFeat.geometry()
            minusBuffer = geom.buffer(
                minusBufferLength, -1, Qgis.EndCapStyle.Round, Qgis.JoinStyle.Round, -1
            )
            if minusBuffer.isEmpty():
                continue
            localHilltopLyr = layerHandler.createMemoryLayerWithFeature(
                hillTopsLyr, hilltopFeat, context
            )
            clippedRaster = self.algRunner.runClipRasterLayer(
                rasterLyr,
                mask=localHilltopLyr,
                context=context,
                nodata=-9999,
                outputRaster=QgsProcessingUtils.generateTempFilename(
                    f"local_clip_{str(uuid4().hex)}.tif"
                ),
            )
            clippedRasterLyr = QgsProcessingUtils.mapLayerFromString(
                clippedRaster, context
            )
            if clippedRasterLyr is None:
                continue
            newFeat = rasterHandler.createMaxPointFeatFromRasterLayer(
                inputRaster=clippedRasterLyr,
                fields=fields,
                fieldName="cota",
                defaultAtributeMap=dict(self.defaultAttrMap),
            )
            featList.append(newFeat)

            if multiStepFeedback is not None:
                multiStepFeedback.setProgress(current * stepSize)
        return featList, clippedHilltops

    def extractAllHilltopPoints(
        self,
        contourLyr: QgsVectorLayer,
        polygonLyr: QgsVectorLayer,
        rasterLyr: QgsRasterLayer,
        geographicBoundsLyr: QgsVectorLayer,
        minusBuffergeographicBoundsLyr: QgsVectorLayer,
        originEpsg: QgsCoordinateReferenceSystem,
        fields: QgsFields,
        context: QgsProcessingContext,
        feedback: QgsFeedback = None,
    ) -> List[QgsFeature]:
        """Extract all hilltop points (both normal and depression) without
        applying distance/grid/exclusion filtering. Unlike
        extractElevationPointsFromHilltops, this includes hilltops of any
        order count (including order 1)."""
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(3, feedback)
            if feedback is not None
            else None
        )
        spatialRelationsHandler = SpatialRelationsHandler()
        layerHandler = LayerHandler()
        currentStep = 0
        if multiStepFeedback is not None:
            multiStepFeedback.setCurrentStep(currentStep)
            multiStepFeedback.pushInfo(self.tr("Finding all hilltops..."))
        hillTopsLyr = spatialRelationsHandler.createHilltopLayerFromPolygonLayer(
            polygonLayer=polygonLyr,
            geographicBoundsLyr=geographicBoundsLyr,
            context=context,
            feedback=multiStepFeedback,
            computeOrder=True,
        )
        self.algRunner.runCreateSpatialIndex(
            hillTopsLyr, context=context, is_child_algorithm=True
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        minusBufferLength = geometryHandler.convertDistance(
            self.contourBufferLength,
            originEpsg=originEpsg,
            destinationEpsg=geographicBoundsLyr.crs(),
        )
        clippedHilltops = self.algRunner.runClip(
            inputLayer=hillTopsLyr,
            overlayLayer=minusBuffergeographicBoundsLyr,
            context=context,
            feedback=multiStepFeedback,
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        nFeats = clippedHilltops.featureCount()
        if nFeats == 0:
            return []
        self.algRunner.runCreateSpatialIndex(
            clippedHilltops, context=context, is_child_algorithm=True
        )
        clippedHilltops = self.algRunner.runJoinAttributesByLocation(
            inputLyr=clippedHilltops,
            joinLyr=contourLyr,
            context=context,
            discardNonMatching=True,
            method=1,
        )
        stepSize = 100 / nFeats
        featList = []
        hilltopList = sorted(
            filter(
                lambda x: x["order_count"] > 1,
                clippedHilltops.getFeatures(),
            ),
            key=lambda x: (x["order_count"], 1.0 / max(x.geometry().area(), 1e-15)),
            reverse=True,
        )
        for current, hilltopFeat in enumerate(hilltopList):
            if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                break
            geom = hilltopFeat.geometry()
            minusBuffer = geom.buffer(
                minusBufferLength,
                -1,
                Qgis.EndCapStyle.Round,
                Qgis.JoinStyle.Round,
                -1,
            )
            if minusBuffer.isEmpty():
                continue
            localHilltopLyr = layerHandler.createMemoryLayerWithFeature(
                hillTopsLyr, hilltopFeat, context
            )
            clippedRaster = self.algRunner.runClipRasterLayer(
                rasterLyr,
                mask=localHilltopLyr,
                context=context,
                nodata=-9999,
                outputRaster=QgsProcessingUtils.generateTempFilename(
                    f"local_clip_{str(uuid4().hex)}.tif"
                ),
            )
            clippedRasterLyr = QgsProcessingUtils.mapLayerFromString(
                clippedRaster, context
            )
            if clippedRasterLyr is None:
                continue
            newFeat = rasterHandler.createMaxPointFeatFromRasterLayer(
                inputRaster=clippedRasterLyr,
                fields=fields,
                fieldName="cota",
                defaultAtributeMap=dict(self.defaultAttrMap),
            )
            featList.append(newFeat)
            if multiStepFeedback is not None:
                multiStepFeedback.setProgress(current * stepSize)
        return featList

    def keepThirdOrderOrHigherPoints(
        self,
        points: Union[List, QgsVectorLayer],
        linesLyr: QgsVectorLayer,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> List[QgsFeature]:
        iterator = (
            points.getFeatures() if isinstance(points, QgsVectorLayer) else points
        )
        nFeats = (
            points.featureCount() if isinstance(points, QgsVectorLayer) else len(points)
        )
        if nFeats == 0:
            return []
        featDict = dict()
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(1, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        stepSize = 100 / nFeats
        for current, feat in enumerate(iterator):
            if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                break
            geom = feat.geometry()
            geomWkb = geom.asWkb()
            if geomWkb in featDict:
                continue
            buffer = geom.buffer(1e-5, -1)
            bbox = buffer.boundingBox()
            nFeats = sum(
                1
                for _ in (
                    i
                    for i in linesLyr.getFeatures(bbox)
                    if i.geometry().intersects(buffer)
                )
            )
            if nFeats < 3:
                continue
            featDict[geomWkb] = feat  # uses dict to remove duplicated features
            if multiStepFeedback is not None:
                multiStepFeedback.setProgress(current * stepSize)

        return list(featDict.values())

    def getElevationPointsFromPlaneAreas(
        self,
        contourLyr: QgsVectorLayer,
        geographicBoundsLyr: QgsVectorLayer,
        rasterLyr: QgsRasterLayer,
        exclusionLyr: QgsVectorLayer,
        polygonLyr: QgsVectorLayer,
        originEpsg: QgsCoordinateReferenceSystem,
        bufferDistance: float,
        contourAreaDict: Dict,
        contourHeightInterval: float,
        gridLyr: QgsVectorLayer,
        gridDict: Dict[QByteArray, int],
        maxPointsPerGridUnit: int,
        fields: QgsFields,
        context: QgsProcessingContext,
        feedback: QgsFeedback,
    ) -> List[QgsFeature]:
        layerHandler = LayerHandler()
        multiStepFeedback = (
            QgsProcessingMultiStepFeedback(4, feedback)
            if feedback is not None
            else None
        )
        if multiStepFeedback is not None:
            currentStep = 0
            multiStepFeedback.setCurrentStep(currentStep)
        localGridSize = geometryHandler.convertDistance(
            self.planeGridSpacing,
            originEpsg=originEpsg,
            destinationEpsg=geographicBoundsLyr.crs(),
        )
        planeGrid = self.algRunner.runCreateGrid(
            extent=geographicBoundsLyr.extent(),
            crs=geographicBoundsLyr.crs(),
            hSpacing=localGridSize,
            vSpacing=localGridSize,
            context=context,
            feedback=multiStepFeedback,
        )
        self.algRunner.runCreateSpatialIndex(
            planeGrid, context, is_child_algorithm=True
        )
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        candidateGridLyr = self.algRunner.runExtractByLocation(
            inputLyr=planeGrid,
            intersectLyr=contourLyr,
            predicate=[2],
            context=context,
            feedback=multiStepFeedback,
        )
        nFeats = candidateGridLyr.featureCount()
        if nFeats == 0:
            return []
        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        stepSize = 100 / nFeats
        pointList = []
        for current, feat in enumerate(candidateGridLyr.getFeatures()):
            if multiStepFeedback is not None and multiStepFeedback.isCanceled():
                break
            localGridLyr = layerHandler.createMemoryLayerWithFeature(
                candidateGridLyr, feat, context
            )
            clippedRaster = self.algRunner.runClipRasterLayer(
                rasterLyr,
                mask=localGridLyr,
                context=context,
                nodata=-9999,
                outputRaster=QgsProcessingUtils.generateTempFilename(
                    f"local_clip_{str(uuid4().hex)}.tif"
                ),
            )
            clippedRasterLyr = QgsProcessingUtils.mapLayerFromString(
                clippedRaster, context
            )
            if clippedRasterLyr is None:
                continue
            newFeatList = rasterHandler.createMaxPointFeatListFromRasterLayer(
                inputRaster=clippedRasterLyr,
                fields=fields,
                fieldName="cota",
            )
            if not newFeatList:
                continue
            tempLayer = layerHandler.createMemoryLayerWithFeatures(
                featList=newFeatList,
                fields=fields,
                crs=rasterLyr.crs(),
                wkbType=QgsWkbTypes.Point,
                context=context,
            )

            # Apply distance-based filtering within this grid cell
            localMinDistance = bufferDistance * 2
            spatiallyFilteredFeats = self.filterFeaturesByDistanceAndExclusionLayer(
                candidatesPointLyr=tempLayer,
                exclusionLyr=exclusionLyr,
                distance=localMinDistance,
                context=context,
                feedback=None,
            )

            pointList.extend(spatiallyFilteredFeats)
            if multiStepFeedback is not None:
                multiStepFeedback.setProgress(current * stepSize)

        if multiStepFeedback is not None:
            currentStep += 1
            multiStepFeedback.setCurrentStep(currentStep)
        return self.filterWithAllCriteria(
            inputPointList=pointList,
            referenceLyr=polygonLyr,
            exclusionLyr=exclusionLyr,
            polygonLyr=polygonLyr,
            bufferDistance=bufferDistance,
            contourAreaDict=contourAreaDict,
            contourHeightInterval=contourHeightInterval,
            gridLyr=gridLyr,
            gridDict=gridDict,
            maxPointsPerGridUnit=maxPointsPerGridUnit,
            fields=fields,
            context=context,
            feedback=multiStepFeedback,
        )

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

    def createInstance(self):
        return ExtractElevationPoints()

    def name(self):
        return "extractelevationpoints"

    def displayName(self):
        return self.tr("Extract Spot Elevation")

    def group(self):
        return self.tr("QA Tools: Terrain Processes")

    def groupId(self):
        return "DSGTools - QA Tools: Terrain Processes"

    def shortHelpString(self):
        return self.tr(
            "This algorithm extracts elevation points from DEM. The output crs is always 4674."
        )
