Estadística zonal con puntos espaciados regularmente dentro de los rasgos (features) de un vectorial tipo polígono con PyQGIS

En el post anterior se refirió la manera de cómo crear puntos espaciados regularmente dentro de los rasgos (features) de un vectorial tipo polígono con PyQGIS. Vamos a aprovechar este hecho para, con un ráster base como capa subyacente, tome los valores de la réticula de éste y, mediante el método ‘identify’ de QgsRasterDataProvider muestree todos los values por feature con el fin de determinar parte de sus estadísticas locales (en este caso número de valores, media y desviación estándar).

El código completo es el siguiente:

import numpy as np

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feats = [ feat for feat in layers[0].getFeatures() ]

#raster layer
xsize = layers[1].rasterUnitsPerPixelX()
ysize = layers[1].rasterUnitsPerPixelY()
provider = layers[1].dataProvider()

points = []
values = [ [] for i in range(len(feats)) ]

print "Wait..."

for k, feat in enumerate(feats):
    extent = feat.geometry().boundingBox()
    xmin = extent.xMinimum()
    ymax = extent.yMaximum()
    xmax = extent.xMaximum()
    ymin = extent.yMinimum()
    
    rows = int((ymax - ymin)/ysize)
    cols = int((xmax - xmin)/xsize)
    
    x = xmin
    y = ymax
    
    geom_feat = feat.geometry()
    
    for i in range(rows+1):
        for j in range(cols+1):
            pt = QgsPoint(x,y)
            tmp_pt = QgsGeometry.fromPoint(pt)
            if tmp_pt.within(geom_feat):
                value = provider.identify(pt, 
                                          QgsRaster.IdentifyFormatValue).results()[1]

                values[k].append(value)
                points.append(tmp_pt.asPoint())
                
            x += xsize
        x = xmin
        y -= ysize

epsg = layers[0].crs().postgisSrid()

#points
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

for value in values:
    print 'count: {} mean: {:.2f} stdev: {:.2f}'.format(np.size(value), np.mean(value), np.std(value))

Cuando se ejecuta el código anterior se obtiene lo siguiente:

grid1

No obstante, se puede prescindir de la creación de la capa vectorial de puntos porque aquí sólo se usó para verificar que el muestreo ráster se efectuó de la manera esperada. Además, los valores de la estadística zonal que se imprimen en la Python Console se pueden incorporar como campos en el vectorial tipo polígono.

Esta entrada fue publicada en PyQGIS, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s