Estadística zonal con inclusión de píxeles recubiertos sólo parcialmente

La estadística zonal se refiere a la determinación de los parámetros estadísticos de un ráster que subyacen a un vectorial. Cuando el tamaño de estos últimos es comparable al tamaño de celda del ráster entonces los efectos de borde del vectorial que tienden a excluir el centroide de las celdas del ráster pueden llegar a ser importantes porque sus values no serán considerados en las estadísticas.

Para ejemplificar este aspecto, en la imagen siguiente se tiene un archivo vectorial que representa el área buffer de 100 m, alrededor de dos puntos, sobre un ráster de 20 filas x 20 columnas con resolución 30 m x 30 m. El objetivo es desarrollar un código que permita ponderar los píxeles del borde de los círculos con base en el recubrimiento parcial.

El código empleado se encuentra a continuación. Al comienzo, determina los cuatro vértices para cada celda del ráster y verifica si sólo uno de ellos cae en el interior del buffer. Si es así, la geometría de la celda (tipo polígono) es interceptada con el buffer y su centroide se emplea para recuperar el value del ráster mediante el método 'identify' de la clase QgsRasterDataProvider de PyQGIS. El índice de solapamiento (indx_inter) basado en las áreas de las zonas intersectadas y las celdas puede ser usado como factor de ponderación una vez recuperados los valores de celda.

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

mapcanvas = iface.mapCanvas() layers = mapcanvas.layers() feat_buffer = layers[1].getFeatures().next() #first feature of buffer layer prov_raster = layers[2].dataProvider() raster_extent = layers[2].extent() xmin, ymin, xmax, ymax = raster_extent.toRectF().getCoords() rows = layers[2].height() cols = layers[2].width() sizeX = layers[2].rasterUnitsPerPixelX() sizeY = layers[2].rasterUnitsPerPixelY() points = [] x = xmin for i in range(rows): for j in range(cols): points.append([QgsPoint(x + i*sizeX, ymax + j*(-sizeY)), QgsPoint(x + (i+1)*sizeX, ymax + j*(-sizeY)), QgsPoint(x + (i+1)*sizeX, ymax + (j+1)*(-sizeY)), QgsPoint(x + i*sizeX, ymax + (j+1)*(-sizeY))]) x = xmin for item in points: geom1 = feat_buffer.geometry() if QgsGeometry.fromPoint(item[0]).within(geom1) or \ QgsGeometry.fromPoint(item[1]).within(geom1) or \ QgsGeometry.fromPoint(item[2]).within(geom1) or \ QgsGeometry.fromPoint(item[3]).within(geom1): item.append(item[0]) geom2 = QgsGeometry.fromPolygon([item]) #cell geom3 = geom2.intersection(geom1) #intersection cell and buffer indx_inter = geom3.area()/geom2.area() centroid = geom3.centroid().asPoint() value = prov_raster.identify(centroid, QgsRaster.IdentifyFormatValue).results()[1] print indx_inter, value

La ejecución del código anterior en la Python Console de QGIS produce los siguientes valores. Si el índice de solapamiento es igual a 1 señala que el recubrimiento de la celda por el shapefile es completo.

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

0.0625204194487 5.0 0.233286907572 5.0 0.128285659439 7.0 0.430554583203 9.0 0.964620626089 7.0 1.0 3.0 0.998752162581 8.0 0.670804885299 6.0 0.0382110109824 6.0 0.259108201581 10.0 0.997982751942 8.0 1.0 2.0 1.0 7.0 1.0 7.0 1.0 6.0 0.569683803093 7.0 0.610928539309 8.0 1.0 2.0 1.0 10.0 1.0 9.0 1.0 2.0 1.0 7.0 0.923521388879 8.0 0.636179042117 9.0 1.0 4.0 1.0 2.0 1.0 2.0 1.0 10.0 1.0 1.0 0.948771891687 10.0 0.341485379809 4.0 1.0 3.0 1.0 8.0 1.0 4.0 1.0 4.0 1.0 8.0 0.654078229379 7.0 0.00328203975005 6.0 0.592744891916 1.0 0.999302262394 3.0 1.0 6.0 1.0 3.0 0.817116132863 4.0 0.0908059107671 3.0 0.192046339914 3.0 0.397494464342 5.0 0.29124537879 6.0 0.0178963092058 9.0

En la imagen siguiente se ejemplifica cómo actua el proceso de iteración empleando sólo el último polígono de intersección y su centroide; que han sido visualizados gracias al plugin QuickWKT. El plugin Value Tool permite corroborar que la selección del value (pixel value igual a 9) fue realizada de la manera esperada para este polígono.