PyQGIS Script – Solving Feature Update Issue in QGIS Script Where Only First Feature is Updated

pyqgis

CONTEXT

I am writing a script to help me to calculate volumes based om DEM for specific areas. The script takes a polygon layer with one or more features, as well as a DEM. There is no output, but a field with calculated volumes for each area should be added.

ISSUE

After finished, the new field is added, but only the first feature has the volume value saved. The others remain NULL

The files I'm using can be downloaded here:

Here's the code:

from qgis.core.additions.edit import edit
from PyQt5.QtCore import (QCoreApplication,
                          QVariant)
from qgis.core import (QgsProcessing,
                       QgsFields,
                       QgsField,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterVectorLayer,
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterField,
                       QgsProcessingParameterEnum,
                       QgsProcessingParameterString,
                       QgsFeatureSink,
                       QgsExpression)
import processing


class ThicknessCalculatorCustomAreas(QgsProcessingAlgorithm):
    POLYGON = 'POLYGON'
    RASTER = 'RASTER'
    CALC_METHOD = 'CALC_METHOD'

    def tr(self, string):

        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return ThicknessCalculatorCustomAreas()

    def name(self):
        return 'Thickness Calculator (Custom areas)'

    def displayName(self):
        return 'Thickness Calculator (Custom areas)'

    def group(self):
        return 'Custom'

    def groupId(self):
        return 'Custom'

    def shortHelpString(self):
        return """""" # TODO

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.POLYGON,
                'Vector file with delimited areas (Only polygon)',
                types=[QgsProcessing.TypeVectorPolygon],
                defaultValue=None))

        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.RASTER,
                'Raster file',
                defaultValue=None))

        self.addParameter(
            QgsProcessingParameterEnum(
                self.CALC_METHOD,
                self.tr('Calculation method:'),
                options=['[0] Count Only Above Base Level: only pixels above the base level will add to the volume.',
                         '[1] Count Only Below Base Level: only pixels below the base level will add to the volume',
                         '[2] Subtract Volumes Below Base level: pixels above the base level will add to the volume, pixels below the base level will subtract from the volume.',
                         '[3] Add Volumes Below Base level: Add the volume regardless whether the pixel is above or below the base level. This is equivalent to sum the absolute values of the difference between the pixel value and the base level.'],
                allowMultiple=False,
                usesStaticStrings=False,
                defaultValue=0))

    def processAlgorithm(self, parameters, context, feedback):

        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        results = {}
        outputs = {}

        polygon = self.parameterAsVectorLayer(
            parameters,
            self.POLYGON,
            context)

        raster = self.parameterAsRasterLayer(
            parameters,
            self.RASTER,
            context)

        calc_method = self.parameterAsEnum(
            parameters,
            self.CALC_METHOD,
            context)

        polygon.dataProvider().addAttributes([QgsField("Calculated_Volume", QVariant.Double)])
        polygon.updateFields()

        fid = polygon.fields()[0]

        with edit(polygon):
            for feature in polygon.getFeatures():
                feat = next(polygon.getFeatures())
                # Extract by attribute
                alg_params = {
                    'FIELD': fid.name(),
                    'INPUT': polygon,
                    'OPERATOR': calc_method,
                    'VALUE': feat[fid.name()],
                    'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
                }
                outputs['ExtractByAttribute'] = processing.run('native:extractbyattribute', alg_params, context=context,
                                                               feedback=feedback, is_child_algorithm=True)
    
                if feedback.isCanceled():
                    return {}
    
                # Clip raster by mask layer
                alg_params = {
                    'ALPHA_BAND': False,
                    'CROP_TO_CUTLINE': True,
                    'DATA_TYPE': 0,  # Use Input Layer Data Type
                    'EXTRA': '',
                    'INPUT': raster,
                    'KEEP_RESOLUTION': False,
                    'MASK': outputs['ExtractByAttribute']['OUTPUT'],
                    'MULTITHREADING': False,
                    'NODATA': None,
                    'OPTIONS': '',
                    'SET_RESOLUTION': False,
                    'SOURCE_CRS': raster,
                    'TARGET_CRS': raster,
                    'X_RESOLUTION': None,
                    'Y_RESOLUTION': None,
                    'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
                }
                outputs['ClipRasterByMaskLayer'] = processing.run('gdal:cliprasterbymasklayer', alg_params, context=context,
                                                                  feedback=feedback, is_child_algorithm=True)
    
                if feedback.isCanceled():
                    return {}
    
                # Raster surface volume
                alg_params = {
                    'BAND': 1,
                    'INPUT': outputs['ClipRasterByMaskLayer']['OUTPUT'],
                    'LEVEL': 0,
                    'METHOD': 0,  # Count Only Above Base Level
                    'OUTPUT_TABLE': QgsProcessing.TEMPORARY_OUTPUT
                }
                outputs['RasterSurfaceVolume'] = processing.run('native:rastersurfacevolume', alg_params, context=context,
                                                                feedback=feedback, is_child_algorithm=True)
    
                if feedback.isCanceled():
                    return {}
    
                feat["Calculated_Volume"] = round(outputs['RasterSurfaceVolume']['VOLUME'], 1)
                polygon.updateFeature(feat)

          return results

QUESTION

What should I change so that all features will be updated with the expected values?

Best Answer

Without testing your entire script, the problem you describe is because feat = next(polygon.getFeatures()) is only storing the first feature, therefore subsequent references to feat are acting on the first feature not the feature loop variable. The simplest fix would be to delete that line and change the loop variable from feature to feat e.g.

for feat in polygon.getFeatures():
    # Extract by attribute

Alternatively, if you want to work with a while loop, you could use the nextFeature() method of QgsFeatureIterator.

E.g.

it = polygon.getFeatures()
feat = QgsFeature()
while it.nextFeature(feat):
    # Do something with feat
Related Question