I would like to extract and convert to a shapefile the coordinates of the highest pixel values located in several polygons. Such treatment was quite straight forward to achieve in R but I'm a total newbie at python and after several unsuccessful attemps I'm asking for some help...

One of the difficulties is that I don't want to use python modules that are not implemented in pyqgis since I wanna remain in the pyqgis console...

EDIT (June 9th)

I digged around to find a solution within Qgis and thanks to this suggestion, I came up with a solution to my question.

Here is the code:

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsProcessingParameterFeatureSink
from qgis.core import QgsProcessingParameterRasterDestination
from qgis.core import QgsProcessingParameterVectorDestination
import processing

class Modele_etape1(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterRasterLayer('mnt', 'MNT', defaultValue=None))
        self.addParameter(QgsProcessingParameterVectorLayer('zoneshumides', 'Zones Humides', types=[QgsProcessing.TypeVector], defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Zh_tampon_50m', 'ZH_tampon_50m', type=QgsProcessing.TypeVectorPolygon, createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Zh_flow_acc', 'ZH_flow_acc', optional=True, type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Exutoires', 'exutoires', type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterRasterDestination('Catchment_area', 'catchment_area', createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterRasterDestination('MntRempli', 'MNT rempli', createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterVectorDestination('Points', 'points', type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(10, model_feedback)
        results = {}
        outputs = {}

        # De polygones à lignes
        alg_params = {
            'INPUT': parameters['zoneshumides'],
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        outputs['DePolygonesLignes'] = processing.run('qgis:polygonstolines', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        if feedback.isCanceled():
            return {}

        # Fill sinks (wang & liu)
        alg_params = {
            'ELEV': parameters['mnt'],
            'MINSLOPE': 0.01,
            'FDIR': QgsProcessing.TEMPORARY_OUTPUT,
            'FILLED': parameters['MntRempli'],
            'WSHED': QgsProcessing.TEMPORARY_OUTPUT
        outputs['FillSinksWangLiu'] = processing.run('saga:fillsinkswangliu', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['MntRempli'] = outputs['FillSinksWangLiu']['FILLED']

        if feedback.isCanceled():
            return {}

        # Tampon
        alg_params = {
            'DISSOLVE': False,
            'DISTANCE': 50,
            'END_CAP_STYLE': 0,
            'INPUT': outputs['DePolygonesLignes']['OUTPUT'],
            'JOIN_STYLE': 0,
            'MITER_LIMIT': 2,
            'SEGMENTS': 5,
            'OUTPUT': parameters['Zh_tampon_50m']
        outputs['Tampon'] = processing.run('native:buffer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Zh_tampon_50m'] = outputs['Tampon']['OUTPUT']

        if feedback.isCanceled():
            return {}

        # Catchment area
        alg_params = {
            'ELEVATION': outputs['FillSinksWangLiu']['FILLED'],
            'METHOD': 0,
            'FLOW': parameters['Catchment_area']
        outputs['CatchmentArea'] = processing.run('saga:catchmentarea', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Catchment_area'] = outputs['CatchmentArea']['FLOW']

        if feedback.isCanceled():
            return {}

        # Zonal statistics
        alg_params = {
            'COLUMN_PREFIX': '_',
            'INPUT_RASTER': outputs['CatchmentArea']['FLOW'],
            'INPUT_VECTOR': outputs['Tampon']['OUTPUT'],
            'RASTER_BAND': 1,
            'STATISTICS': [6]
        outputs['ZonalStatistics'] = processing.run('native:zonalstatistics', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        if feedback.isCanceled():
            return {}

        # Découper un raster selon une couche de masque
        alg_params = {
            'ALPHA_BAND': False,
            'CROP_TO_CUTLINE': True,
            'DATA_TYPE': 0,
            'EXTRA': '',
            'INPUT': outputs['CatchmentArea']['FLOW'],
            'KEEP_RESOLUTION': False,
            'MASK': outputs['Tampon']['OUTPUT'],
            'MULTITHREADING': False,
            'NODATA': None,
            'OPTIONS': '',
            'SET_RESOLUTION': False,
            'SOURCE_CRS': 'ProjectCrs',
            'TARGET_CRS': 'ProjectCrs',
            'X_RESOLUTION': None,
            'Y_RESOLUTION': None,
            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
        outputs['DcouperUnRasterSelonUneCoucheDeMasque'] = processing.run('gdal:cliprasterbymasklayer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        if feedback.isCanceled():
            return {}

        # Raster values to points
        alg_params = {
            'GRIDS': outputs['CatchmentArea']['FLOW'],
            'NODATA': True,
            'POLYGONS': outputs['ZonalStatistics']['INPUT_VECTOR'],
            'TYPE': 0,
            'SHAPES': parameters['Points']
        outputs['RasterValuesToPoints'] = processing.run('saga:rastervaluestopoints', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Points'] = outputs['RasterValuesToPoints']['SHAPES']

        if feedback.isCanceled():
            return {}

        # Join attributes by location
        alg_params = {
            'DISCARD_NONMATCHING': False,
            'INPUT': outputs['RasterValuesToPoints']['SHAPES'],
            'JOIN': outputs['Tampon']['OUTPUT'],
            'JOIN_FIELDS': ['SITE_CODE'],
            'METHOD': 0,
            'PREDICATE': [0],
            'PREFIX': '',
            'OUTPUT': parameters['Zh_flow_acc']
        outputs['JoinAttributesByLocation'] = processing.run('native:joinattributesbylocation', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Zh_flow_acc'] = outputs['JoinAttributesByLocation']['OUTPUT']

        if feedback.isCanceled():
            return {}

        # Sélection par expression
        alg_params = {
            'EXPRESSION': '\"FLOW\" = maximum(\"FLOW\", \"SITE_CODE\")',
            'INPUT': outputs['JoinAttributesByLocation']['OUTPUT'],
            'METHOD': 0
        outputs['SlectionParExpression'] = processing.run('qgis:selectbyexpression', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        if feedback.isCanceled():
            return {}

        # Extraire les entités sélectionnées
        alg_params = {
            'INPUT': outputs['SlectionParExpression']['OUTPUT'],
            'OUTPUT': parameters['Exutoires']
        outputs['ExtraireLesEntitsSlectionnes'] = processing.run('native:saveselectedfeatures', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Exutoires'] = outputs['ExtraireLesEntitsSlectionnes']['OUTPUT']
        return results

    def name(self):
        return 'Modele_etape1'

    def displayName(self):
        return 'Modele_etape1'

    def group(self):
        return ''

    def groupId(self):
        return ''

    def createInstance(self):
        return Modele_etape1()

The only remaining problem is that I can only run this script / model in QGIS 3.12...

I tried to use it in QGIS 3.10 but the tool "join by location" does not work within my model (it works if use it as a single tool in a QGIS project though...)

Does someone know if it is a problem related to this specific Qgis version ? Is there a way to fix it ?

Many thanks in advance !

  • You passed a variable of type str to QgsZonalStatistics. Do you know what type it was expecting? – Paul H Jun 5 at 15:15
  • Hey Paul ! I don't understand why I get such error message since I'm quite sure I have entered the right variables to QgsZonalStatistics qgis.org/api/classQgsZonalStatistics.html – wanderzen Jun 5 at 15:20
  • Based on that documentation, I would expect the method to require some sort of Layer object, not a string representation of the file path of the data – Paul H Jun 5 at 16:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.