Manejando píxeles con solapamiento parcial para realizar estadística zonal con PyQGIS en QGIS 3

El manejo de píxeles con solapamiento parcial para realizar estadística zonal robusta con PyQGIS ya había sido tratado en esta pregunta de https://gis.stackexchange.com. No obstante, el planteamiento en ésta otra es algo diferente pero el código puede adaptarse para satisfacer el requerimiento. 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 54 55 56 57 58 59 60 61 62 63 64 65 66 67

import processing registry = QgsProject.instance() Polygons = registry.mapLayersByName('Polygons_new') iterator = Polygons[0].getFeatures() feats_Polygons = [ feat for feat in iterator ] Raster = registry.mapLayersByName('Raster') rprovider = Raster[0].dataProvider() xmin, ymin, xmax, ymax = Raster[0].extent().toRectF().getCoords() extent = str(xmin)+ ',' + str(xmax)+ ',' +str(ymin)+ ',' +str(ymax) parameters = {"TYPE":2, "EXTENT":extent, "HSPACING":1000, "VSPACING":1000, "HOVERLAY":0, "VOVERLAY":0, "CRS":Raster[0].crs().authid(), "OUTPUT":"/home/zeito/pyqgis_data/tmp_grid.shp"} path = processing.run('qgis:creategrid', parameters) grid = QgsVectorLayer(path['OUTPUT'], 'grid', 'ogr') feats_grid = [ feat for feat in grid.getFeatures() ] geom_pol = [] attr = [] for featg in feats_grid: for featp in feats_Polygons: if featg.geometry().intersects(featp.geometry()): geom = featg.geometry().intersection(featp.geometry()) geom_pol.append(geom.asWkt()) pt = geom.centroid().asPoint() value = rprovider.identify(pt, QgsRaster.IdentifyFormatValue).results()[1] attr.append([value, geom.area()]) epsg = Polygons[0].crs().postgisSrid() uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=value:double&field=area:double&field=w_value:double""&index=yes" mem_layer = QgsVectorLayer(uri, 'polygon', 'memory') prov = mem_layer.dataProvider() feats = [ QgsFeature() for i in range(len(geom_pol)) ] for i, feat in enumerate(feats): feat.setAttributes([i, attr[i][0], attr[i][1],(attr[i][0]*(attr[i][1]/1e6))]) feat.setGeometry(QgsGeometry.fromWkt(geom_pol[i])) prov.addFeatures(feats) QgsProject.instance().addMapLayer(mem_layer)

Ejecutado en la Python Console de QGIS 3, para el shapefile de la imagen siguiente:

permite obtener el resultado que se presenta a continuación:

Se puede observar que el solapamiento es proporcional al área barrida por el vectorial en la región que solapa cada píxel.

En el ejemplo siguiente es aún más fácil la corroboración porque existe un píxel totalmente cubierto (100 % de solapamiento) y, por tanto, el factor de ponderación es 1 y, en consecuencia, el píxel value coincide con el valor ponderado.