Estadística zonal robusta con PyQGIS

Uno de los problemas con las herramientas de estadística zonal que vienen generalmente incorporadas con los SIGs es que producen resultados aproximados cuando los rasgos (features) son mucho menores que las celdas del ráster. Esto se debe a que se pueden solapar hasta 4 de ellas, con diferente grado, y no lo ponderan sino que tienden a considerar el de mayor proporción ignorando los otros.

Esto hace que la estadística no sea robusta y se procuren métodos con ejecución de código personalizado para solventar este problema.

El código a continuación ha sido desarrollado para llevar a cabo la estadística zonal del solapamiento de los rasgos 19 a 27 del vectorial y el ráster que se encuentran en este link. La imagen con zoom a los rasgos se coloca a continuación.

La visión general de todos los rasgos del shapefile lucen como a continuación:

Cuando se ejecuta el código siguiente en la Python Console de QGIS:

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

import processing registry = QgsMapLayerRegistry.instance() Polygons = registry.mapLayersByName('Polygons') request = QgsFeatureRequest().setFilterExpression (u'"ID_sub.pat" >= 19 AND "ID_sub.pat" <= 27') iterator = Polygons[0].getFeatures(request) 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":1, "EXTENT":extent, "HSPACING":1000, "VSPACING":1000, "CRS":Raster[0].crs().authid(), "OUTPUT":None} path = processing.runalg('qgis:creategrid', parameters) grid = QgsVectorLayer(path['OUTPUT'], 'grid', 'ogr') feats_grid = [ feat for feat in grid.getFeatures() ] for featg in feats_grid: for featp in feats_Polygons: if featg.geometry().intersects(featp.geometry()): geom = featg.geometry().intersection(featp.geometry()) pt = geom.centroid().asPoint() value = rprovider.identify(pt, QgsRaster.IdentifyFormatValue).results()[1] print value, geom.area()/62500.0, featp.id()+1

se obtiene este resultado:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

-88.0 0.11367261864 22 -88.0 0.744772910144 25 -86.0 0.0389546011361 22 -86.0 0.152627219776 23 -86.0 0.152627219776 24 -86.0 0.255227089856 25 -86.0 1.0 26 -86.0 1.0 27 -102.399993896 0.744772910144 19 -102.399993896 0.631100291505 22 -92.799987793 0.255227089856 19 -92.799987793 1.0 20 -92.799987793 1.0 21 -92.799987793 0.21627248872 22 -92.799987793 0.847372780224 23 -92.799987793 0.847372780224 24

donde en la primera columna se tienen los raster value, en la segunda columna el factor de ponderación por superficie y en la tercera columna el id del rasgo. Se puede observar que para los polígonos 20,21,26,27 (primera imagen) sus factores de ponderación son uno (tal como se espera). En los otros casos es proporcional a la intersección de los rasgos de los polígonos con la rejilla que representa el ráster.

Los resultados fueron corroborados con Value Tool y QuickWKT QGIS plugins y fueron tal como se esperaba.