Compute zonal statistics in QGIS with overlapping features

Question

I would like to compute zonal statistics in QGIS (preferably using the Python console) of a polygon layer with overlapping features.

I've seen different questions (and answers) on how to compute zonal statistics, but none of them refer specifically to how to deal with overlapping features in QGIS. However, I've read in different posts that overlapping features might complicate things when using, for example, ArcGIS.

The question is:

Do you know if the QgsZonalStatistics class works well with layers with overlapping polygons?

If not, is there any alternative to carry out this task in QGIS? For example, to compute zonal statistics iterating over the different features of the layer?

Answer 1

For computing zonal statistics in PyQGIS with overlapping features you can use 'unary_union' and 'polygonize' methods from shapely python module. For trying out my approach I used next shapefile (with overlapping features) and raster:

Code used is:

import fiona
from shapely.geometry import shape, LineString
from shapely.ops import unary_union, polygonize
from qgis.analysis import QgsZonalStatistics

registry = QgsMapLayerRegistry.instance()

overlapped_features = registry.mapLayersByName('overlapped_features')
raster = registry.mapLayersByName('utah_demUTM2')

vector_provider = overlapped_features[0].dataProvider()

path = vector_provider.dataSourceUri()

path = path.split('|')[0]

c = fiona.open(path)

collection = [ shape(item['geometry']) for item in c ]

rings = [ LineString(pol.exterior.coords) for pol in collection ]

union = unary_union(rings)

new_intersections = [ geom.wkt for geom in polygonize(union) ]

epsg = overlapped_features[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'new_polygon',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(new_intersections)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(new_intersections[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

raster_provider = raster[0].dataProvider()

path_raster = raster_provider.dataSourceUri()

zoneStat = QgsZonalStatistics(mem_layer, path_raster,"", 1, QgsZonalStatistics.All)

zoneStat.calculateStatistics(None)

After running it at the Python Console of QGIS I got:

Resulting layer ('new_polygon') has individualized each overlapping areas and at the image was selected one of them. It can be observed that in all cases, there was computed all possible statistics of QgsZonalStatistics class.

Editing Note:

Based in Vincent's commentary, the answer is yes. It works well for each individual polygon in the original 'overlapped_features' layer, even when there are overlapped features. In this case the code is simpler:

from qgis.analysis import QgsZonalStatistics

registry = QgsMapLayerRegistry.instance()

overlapped_features = registry.mapLayersByName('overlapped_features')
raster = registry.mapLayersByName('utah_demUTM2')

raster_provider = raster[0].dataProvider()

path_raster = raster_provider.dataSourceUri()

zoneStat = QgsZonalStatistics(overlapped_features[0], path_raster,"", 1, QgsZonalStatistics.All)

zoneStat.calculateStatistics(None)

and after running it at the Python Console of QGIS we have: