¿Obtener mínimo, máximo, media, promedio con gdal_rasterize?

Question

No estoy 100% seguro pero parece que gdal_rasterizar toma por defecto el valor de la primera característica encontrada. ¿Hay una manera de decir que me gustaría obtener mínimo / máximo / promedio / media (/ recuento / suma sólo para conseguir realmente de lujo) de las características del vector valores dentro de la celda?

Particularmente tengo un archivo de puntos csv/shapefile con columnas X,Y,Z y puntos donde a menudo dos o más puntos caen en un pixel. ¿Cómo puedo rasterizar el valor máximo de esos puntos (o mínimo, promedio, media)?

Existen sql "-where" y "-sql select_statement", pero no estoy seguro de si podrían utilizarse en este caso.

Answer 1

Puede utilizar GDAL Python para rasterizar, sin embargo, antes es necesario un poco de Python (en mi caso PyQGIS). Probé mi enfoque creando un pequeño raster de 3x3 (para facilitar la validación de los resultados) y, después, utilicé la herramienta de procesamiento 'Create grid' para obtener un shapefile Grid perfectamente alineado con la capa raster. También se creó otro shapefile, con una cantidad variable de puntos por celdas raster. Este shapefile de puntos tiene un campo (llamado 'field'; ver método 'attribute' en el código de abajo) con valores flotantes aleatorios. Todo el conjunto se puede observar en la siguiente imagen:

El siguiente código selecciona los valores en el campo 'field' del shapefile de puntos y, si están dentro de cada característica del shapefile Grid, se agrupan por índices de fila y columna (para trabajar en el caso de puntos aleatorios sin id consecutivo por celda). Estos valores también se imprimen en la consola Python de QGIS.

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]) )

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

Después de ejecutar el código en la consola Python de QGIS los valores impresos fueron:

[[58.1, 62.66, 12.93], [40.23], [93.21, 41.74], [3.23], [63.45, 56.86, 63.22, 57.62], [22.0], [48.23, 95.04], [50.23, 35.13], [43.48]]

Hay nueve listas porque el raster es de 3 x 3. Cada lista interna tiene valores de puntos por cada celda y su cantidad es igual al número de puntos. Ahora, es fácil extraer el mínimo, máximo o media por celda y rasterizar estos valores utilizando el módulo python GDAL (después de remodelarlos y convertirlos en array con numpy).

Answer 2

Problema similar en el que tenía cero a muchos puntos dentro de cada celda de la cuadrícula y quiero que la suma de un atributo para todas las características. He utilizado gdal_rasterize en la CLI, (gdal v2.1.2). Las opciones estándar resultan en una banda de todos los ceros.

La clave era poner a cero la constante de inicialización y luego pasar la bandera -add para que se acumularan los valores de cada característica.

gdal_rasterize -te <xmin ymin xmax ymax> -tr <xres yres> -of Float32 -ot GTiff -a <targ_attrib> -init 0.0 -add -l <src_layer> src_layer.shp out_grid.tif