Máximos y mínimos locales con el método ‘ReadAsArray’ de GDAL/Python

El método ‘ReadAsArray’ de GDAL/Python permite extraer los valores de los píxeles de un ráster como un bloque para realizar alguna operación útil con ellos (estadísticas, modificación, álgebra de rásteres, clipping, etc). Sin embargo, también se puede recorrer el ráster por bloques con el fin de calcular, por ejemplo, ciertos parámetros estadísticos en ellos. Este sería el caso de la determinación de sus máximos y mínimos locales o por bloque.

El código siguiente va a ser empleado para determinar los máximos y mínimos locales en un ráster de 20 x 20 en bloques de 10 x 10. Por tanto, se tendrán 4 matrices o arrays de valores para los cuales ésto será posible. No obstante, modificando el valor de np se pueden obtener estadísticas para otros bloques cuadrados. Si la dimensión escogida para el bloque no es múltiplo de 2 o de 5 habrá algunos bloques de menor dimensión que no serán contabilizados.

from osgeo import gdal, osr
import numpy as np

driver = gdal.GetDriverByName('GTiff')
filename = "/home/zeito/pyqgis_data/aleatorio.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

cols = dataset.RasterXSize
rows = dataset.RasterYSize

print cols, rows

nb = 10

for i in range(0, rows, nb):
    for j in range(0, cols, nb):
        print "verts:", i, j
        data = band.ReadAsArray(j, i, nb, nb)
        print data
        min = np.min(data)
        print "min: ", min
        max = np.max(data)
        print "max: ", max
        print "maximum"

        try:
            for m, item in enumerate(data):
                for n, element in enumerate(item):
                    if element == max:
                        print m + i, n + j

            print "minimum"

            for m, item in enumerate(data):
                for n, element in enumerate(item):
                    if element == min:
                        print m + i, n + j
        
        except TypeError:
            print "there is no matrix"

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

minimos1

La última matriz que se imprime es la que comienza en el vértice (10,10) y termina en (19,19). Su valor mínimo es 1 y el máximo 100. Lo que se imprime como máximos y mínimos son los índices (fila, columna) correspondientes.

En la imagen siguiente se observa que el código funciona tal como se espera.

minimos2

Con los índices de (fila, columna) es fácil conformar, por ejemplo, memory layers separadamente para los máximos y mínimos locales.

Esta entrada fue publicada en Código Python, GDAL, 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