¿Cómo QGIS del Zonal Estadísticas de la manija que se superponen parcialmente de píxeles?

Question

Hay un par de Stackexchange puestos de hacer esta pregunta, pero ninguno de ellos eran lo suficientemente contestadas: 1, 2.

Yo soy de la escritura de secuencias de comandos de Python que realizan esencialmente las mismas operaciones como Zonal Estadísticas, y he estado usando Zonal Estadísticas para comprobar mi secuencias de comandos de salidas.

He notado que cuando hay muchas células en el polígono que se está consultando, mis resultados coinciden perfectamente, pero hay discrepancias, cuando sólo hay un par de celdas ráster en el polígono. El siguiente es un ejemplo de esto. Yo soy de computación "de la cuenta" para cada polígono.

Prueba 1: Devuelve el recuento es de 0,45, lo que parece más o menos correcto si los píxeles pueden ser tratados marginalmente en lugar de simplemente in/out.

Prueba 2: Resultado Similar; count es 0.39.

Prueba 3: Con una mayor poylgon, ahora vemos que el recuento es exactamente 2. Tenga en cuenta que hay exactamente 2 píxeles de los centros que figuran en el polígono.

Gráfico 4: Aquí, sólo un píxel centro está dentro del polígono. El recuento es 1.008, aunque. Si se excluye el área de polígonos a partir de los píxeles cuyo centro está en el polígono, el resto del área es, obviamente, mayor de 0.008 de un píxel.

A partir de estas exposiciones, parece que Zonales Estadísticas trata estos casos especialmente. Parece a mí que si hay menos de 2 píxeles de centros dentro del polígono, se realiza algún tipo de promedio como polygon_area/raster_cell_area para devolver el recuento. Debo señalar que he 48 de estos polígonos, y todo volvió a cuenta de mayor que o igual a 2 son enteros; fracciones de cuenta sólo se devuelven para <2.

La pregunta es, entonces, ¿cómo son estas fracciones de píxel cuenta incorporado en el otro resultado de las estadísticas por Zonas Estadísticas? Esta pregunta indica un comportamiento extraño con max/min cuando se utilizan grandes cantidades celdas ráster, y de mis propias pruebas devuelve los valores de la media no es congruente con un simple "pixel centro en o fuera del polígono." Para reiterar, me sale exacta acuerdo cuando el subyacente de trama contiene el número de píxeles dentro del polígono. No podía encontrar ninguna documentación acerca de esto, y puede ser importante cuando las personas que nunca han tenido problemas con la Zonal Stats pero el uso de una gruesa trama y siendo conscientes de que existe un comportamiento diferente cuando sólo hay un par de píxeles dentro del polígono.

Answer 1

Puedo manejar parcialmente superpuestas píxeles adaptación de código (para PyQGIS en QGIS 3) en la respuesta a esta pregunta: Robusto Zonal estadísticas en R y QGIS.

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)

Polygons_new referido en el código anterior se puede observar en la siguiente imagen:

Después de ejecutarlo en la Consola de Python de QGIS 3, el resultado fue:

Tabla de atributos se ha ponderado de los valores de píxel para cada interceptó área del polígono (con una rejilla temporal que coincide con celdas ráster).

En el siguiente ejemplo se puede observar un shapefile (Polygons_new2) que se cruza con un completo píxel. Su valor ponderado en la tabla de atributos es el mismo que el valor del píxel (como se esperaba). Resto de los valores de los píxeles son proporcionales a la superficie cubierta por shapefile.