Iterating over the features of a vector layor then performing zonal statistics algorithm on each feature in PyQGIS

Question

I composed this script in PyQGIS to perform the zonal statistics on a vector layer as for all its features at once, it works perfectly, when the script is finished, I inspected the attribute table and indeed it is fine for every feature in the vector layor.

from qgis.analysis import QgsZonalStatistics
from qgis.core import QgsRasterLayer, QgsVectorLayer

v_layer = QgsVectorLayer("gadm40_DEU_4.shp", "gadm", "ogr")

r_layer = QgsRasterLayer("Deutschland.vrt", "dem")


zone_stat = QgsZonalStatistics (v_layer, r_layer,
    attributePrefix = 'dem_',
    rasterBand = 1,
    stats = QgsZonalStatistics.Statistics(
        QgsZonalStatistics.Count|
        QgsZonalStatistics.Max|
        QgsZonalStatistics.Mean|
        QgsZonalStatistics.Median|
        QgsZonalStatistics.Min|
        QgsZonalStatistics.StDev
    )
)

zone_stat.calculateStatistics(None)

My goal is to do the same but not for whole vector layor at once, so iterate over features then do the zonal statistics because after each iteration I would to do something from the result on this feature.

I tried somwthing like this but it didnt work, is there a method to do it, so without writing the result to the vector layor just only return them with each iteration and I can do whatever I want with them?

My try:

from qgis.analysis import QgsZonalStatistics
from qgis.core import QgsRasterLayer, QgsVectorLayer

v_layer = QgsVectorLayer("gadm40_DEU_4.shp", "gadm", "ogr")

r_layer = QgsRasterLayer("Deutschland.vrt", "dem")


for i, feature in enumerate(v_layer.getFeatures()):
    feature_geometry = feature.geometry()
    #this approach doesnt work
    zone_stat = QgsZonalStatistics (feature_geometry, r_layer,
        attributePrefix = 'dem_',
        rasterBand = 1,
        stats = QgsZonalStatistics.Statistics(
            QgsZonalStatistics.Count|
            QgsZonalStatistics.Max|
            QgsZonalStatistics.Mean|
            QgsZonalStatistics.Median|
            QgsZonalStatistics.Min|
            QgsZonalStatistics.StDev
            )
        )
    
    zone_stat.calculateStatistics(None)
    

    feature_name = feature['NAME_4']
    print(feature_name)
    
    break #breaking for test purposes

Answer 1

You can make use of the QgsFeatureSource.materialize() method.

from qgis.analysis import QgsZonalStatistics
from qgis.core import QgsRasterLayer, QgsVectorLayer

v_layer = QgsVectorLayer("gadm40_DEU_4.shp", "gadm", "ogr")
r_layer = QgsRasterLayer("Deutschland.vrt", "dem")

# for convenience, _stats is defined in the order that the Zonal Statistics tool creates the fields
_stats = ['count', 'mean', 'median', 'stdev', 'min', 'max']

# add the neccessary fields to the original layer
with edit(v_layer):
    for s in _stats:
        variant = QVariant.Int if 'count' in s else QVariant.Double
        v_layer.addAttribute(QgsField('dem_' + s, variant))
    
v_layer.updateFields()

# get the field indexes of the newly added fields
fields = v_layer.fields()
field_idxs = [fields.indexFromName('dem_' + s) for s in _stats]

# reference to the data provider of v_layer 
prov = v_layer.dataProvider()

for feature in v_layer.getFeatures():
    # current feature id
    fid = feature.id()
    
    # make a feature request using the current fid
    req = QgsFeatureRequest()
    req.setFilterFid(fid)
    
    # materialize only the filtered feature into a new temporary layer
    filtered_layer = v_layer.materialize(req)
    
    # calculate the zonal statistics on the temporary layer
    zone_stat = QgsZonalStatistics (filtered_layer, r_layer,
        attributePrefix = '_',
        rasterBand = 1,
        stats = QgsZonalStatistics.Statistics(
            QgsZonalStatistics.Count|
            QgsZonalStatistics.Max|
            QgsZonalStatistics.Mean|
            QgsZonalStatistics.Median|
            QgsZonalStatistics.Min|
            QgsZonalStatistics.StDev
            )
        )
        
    zone_stat.calculateStatistics(None)
    
    # get the one feature from the temporary layer
    feat = next(filtered_layer.getFeatures())
    
    # retrieve the calculated statistics
    # the length of the _stats list is added here because the statistics fields are appended onto the existing fields (those created in the loop above)
    calculated_stats = [feat[idx + len(_stats)] for idx in field_idxs]
    
    # create a dictionary of field indexes and calculated statistics
    attributes_to_change = dict(zip(field_idxs, calculated_stats))

    # create an attribute map with the fid as the key
    attribute_map = {fid: attributes_to_change}

    # update the attribute values using the attribute map
    with edit(v_layer):
       updated = prov.changeAttributeValues(attribute_map)
    
    break #breaking for test purposes

Note: I cannot explain why the field order is different. According to to v_layer.fields().names() they are in the same order in the data provider, but not when displayed in the attribute table.