Estadística zonal con rasgos que se solapan usando PyQGIS

Para determinar la estadística zonal en un vectorial cuyos rasgos se solapan sobre un ráster, se pueden usar los métodos 'unary_union' y 'polygonize' del módulo 'shapely' de Python. El uso combinado de estos métodos permite individualizar las zonas solapadas y facilitar el cálculo de las estadísticas. Para probar esta aproximación se usaron el shapefile y el ráster que se viasualizan en la imagen siguiente:

El código empleado se presenta a continuación:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

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)

Después de ejecutado en la Python Console de QGIS la memory layer producida presenta en su tabla atributiva las estadísticas deseadas; tal como se dispone en la imagen siguiente. Allí se ha resaltado una de las zonas de solape que ha sido individualizada durante el proceso.