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 tofeat
are acting on the first feature not thefeature
loop variable. The simplest fix would be to delete that line and change the loop variable fromfeature
tofeat
e.g.Alternatively, if you want to work with a while loop, you could use the nextFeature() method of
QgsFeatureIterator
.E.g.