Rasterizando por valores mínimos, máximos o promedio de un campo en un vectorial de puntos

La situación que se tiene en la imagen siguiente es la de un pequeño ráster de 3 x 3 (resolución 73.9887 m x 73.9887 m) en cuya extensión existen 17 puntos generados aleatoriamente y donde en algunas celdas se observa que hay más de un punto. Además, en la tabla atributiva del vectorial de puntos se observa que al seleccionar los dos primeros, estos aparecen (en amarillo) en celdas no contiguas en el Map Canvas reflejando que no existe ninguna sistematización en la asignación de los id.

rasterize1

Una forma de rasterizar estos puntos por los valores máximos, mínimos o medios de lo que alberga el campo ‘field’ del vectorial de puntos sería seleccionando y agrupando todos aquellos que se encuentren dentro de una misma celda. Esto puede hacerse empleando el método ‘within’ de QgsGeometry si se crea un vectorial tipo rejilla (polígono grid) que coincida y contemple la resolución del ráster. La herramienta para hacerlo de manera expedita existe en Processing (Create Grid); tal como se presenta en la imagen siguiente:

rasterize2

Con los vectoriales de puntos y rejilla puede implementarse entonces un código que registre cada inclusión de puntos por celda y almacene en una lista los valores por los que se desea rasterizar. Como los ids de los puntos no son secuenciales entonces hay que incorporar a los valores extraidos el índice de fila y columna de la celda para que luego puedan ser desplegados en el orden correcto.

El código completo es el siguiente:

import operator

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

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

grid = [ feat for feat in layers[1].getFeatures() ]

rextent = layers[2].extent()
rprovider = layers[2].dataProvider()

xmin = rextent.xMinimum()
ymax = rextent.yMaximum()
xsize = layers[2].rasterUnitsPerPixelX()
ysize = layers[2].rasterUnitsPerPixelY()

elements = []

for i, point in enumerate(points):
    posX = point.geometry().asPoint()[0]
    posY = point.geometry().asPoint()[1]

    row = int((ymax - posY)/ysize)
    col = int((posX - xmin)/xsize)

    for j, item in enumerate(grid):
        if point.geometry().within(item.geometry()):
            elements.append([row, col, point.attribute('field')])

#elements.sort( key=lambda x: (x[0], x[1]) )
elements.sort( key=operator.itemgetter(0, 1) )

n = len(elements)
    
values = [ [] for i in range(len(grid)) ]
    
k = 0
    
values[0].append(elements[0][2])
    
for i in range(n-1):
    if elements[i][0] == elements[i+1][0] and elements[i][1] == elements[i+1][1]:
        values[k].append(elements[i + 1][2])
    else:
        k += 1
        values[k].append(elements[i+1][2])

print values

Cuando se ejecuta en la python Console de QGIS se obtiene:

[[4.43, 59.19], [21.77, 91.34, 26.32], [23.93], 
[21.83, 3.59], [42.65, 17.05], [13.95, 18.35, 35.3], 
[68.37], [49.61, 5.57], [64.51]]

Puede observarse que es una lista de listas con 9 elementos, es decir, equivalente a 3 x 3. Por separado, las listas internas albergan los valores del campo ‘field’ correspondiente en número de: 2, 3, 1, 2, 2, 3, 1, 2, 1; respectivamente. Para comprobar que el código funciona de forma aceptable se seleccionó arbitrariamente, como ejemplo, la segunda lista (en realidad se probó con todas). Por tanto, en la tabla atributiva se marcaron los valores correspondientes y en el Map Canvas se tiene el resaltado de los puntos con ids 0, 4, 16; tal como se presenta en la imagen siguiente:

rasterize3

Por tanto, si se quiere rasterizar por valores máximos, mínimos o medios de los que albergan individualmente las listas internas, es sólo cuestión de seleccionarlos mediante una extensión muy sencilla del código, convertirlos en array con reshape a 3 x 3 (con numpy) antes de grabarlos como ráster con el módulo GDAL de Python.

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