Compute zonal statistics in QGIS with overlapping features

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?

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]))

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:

• Thanks a lot @xunilk !! This is really helpful, and it might become super handy at some point of the task. However, my question was simpler (I'm sorry if I was not clear or precise enough). I wanted to know if QgsZonalStatistics works well for each individual polygon in the original 'overlapped_features' layer, even when there are overlapped features. Feb 25, 2017 at 10:45
• Yes, it works well for each individual polygon in the original 'overlapped_features' layer, even when there are overlapped features. Please, see my Editing Note. Feb 25, 2017 at 11:28