2

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

1 Answer 1

3

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.

enter image description here

7
  • I tried it, nothing written to original vector layer as expected but where can I get the attributes of this filtered_layer? I tried sthg like filtered_layer.attributes() or filtered_layer.AttributeTable but not worked. How can get those attributes of this layer plus the new added attributes by zonal statistics? For example print them simply
    – Khaled
    Feb 22 at 13:07
  • Apologies, I didn't include the part to write the results back to the original layer, only to perform the zonal statistics on an individual feature. I will update my answer.
    – Matt
    Feb 22 at 13:09
  • Having said that, must the zonal statistics be done for one feature at a time? Would it not be easier to calculate it for the whole layer, and then iterate through the features for further processing?
    – Matt
    Feb 22 at 13:14
  • No, because sometimes I have very large vector and raster layer like for a country like Russia, RAM would crash or it hangs. I did this to get attributes, it worked: ... zone_stat.calculateStatistics(None) for feature_in_filtered_layer in filtered_layer.getFeatures(): print(feature_in_filtered_layer.attributes()) break #because it should be one feature in this filtered layer Thank for you, that helped but is it your method memory efficient? Is there no sthg like clear or close to clear memory after each iteration?
    – Khaled
    Feb 22 at 13:23
  • thank you much, your answer helped me even before the new edit. And yes, I noticed the reverse order of the data got from stats and I took that into consideration.
    – Khaled
    Feb 22 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.