3

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?

3

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:

enter image description here

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:

enter image description here

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:

enter image description here

  • 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. – Vincent Feb 25 '17 at 10:45
  • 2
    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. – xunilk Feb 25 '17 at 11:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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