Subregiones ráster usando GDAL en ambiente de PyQGIS

En artículos anteriores (1, 2, 3, 4, 5, 6, 7) he considerado la manera más eficiente para leer un ráster utilizando una scanline (línea de barrido) y la función ‘unpack’ a través de los módulos GDAL-struct de Python. Se ha aplicado especialmente en la obtención de las estadísticas del ráster; ya bien sea sobre la totalidad de sus valores o sobre las filas o columnas. En este post vamos a considerar ahora su uso en el muestreo de subregiones, por lo que la scanline pasará más bien a tener la connotación de un arreglo matricial de valores.

Para poder verificar más fácilmente que el script de Python funciona de la manera esperada preparé un ráster de 29 filas por 29 columnas y asumí, arbitrariamente, que quería muestrear la subregión 3×3 que comienza en la posición (xoff, yoff) = (5,7). Los values del ráster que se encuentran en la imagen inferior fueron corroborados manualmente.

gdal1

El código siguiente permite seleccionar la subregión delimitada anteriormente y determinar algunos parámetros estadísticos básicos (media, mediana, mínimo, máximo) de los valores de la matriz. Toma el paso a la capa activa en la Map View de QGIS a través del provider. La función ‘ReadRaster’ en la scanline permite tomar el arreglo 3×3, con un buffer idéntico de 3×3, a partir del punto (5,7). La amplitud del buffer (3×3 =9) es la que dicta que valores serán desempaquetados finalmente a través la función miembro de struct (en este caso todos los 9 valores).

from osgeo import gdal
import struct
import numpy as np
layer = iface.activeLayer()
provider = layer.dataProvider()
fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}
path= provider.dataSourceUri()
dataset = gdal.Open(path)
band = dataset.GetRasterBand(1)

BandType = gdal.GetDataTypeName(band.DataType)

scanline = band.ReadRaster(5, 7, 3, 3, 3, 3, band.DataType)
values = struct.unpack(fmttypes[BandType] * 9, scanline)

print values
print "mean = ", np.mean(values)
print "median = ", np.median(values)
print "min = ", np.min(values)
print "max = ", np.max(values)

dataset = None

La ejecución del script anterior a través del Editor de la Python Console permite verificar que la matriz considerada en la imagen anterior es la que fue efectivamente muestreada:

Python Console 
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(u'C:/Users/zeito/pyqgis_scripts/scanline.py'.encode('mbcs'))
(1669, 1666, 1661, 1670, 1667, 1665, 1672, 1671, 1669)
mean =  1667.77777778
median =  1669.0
min =  1661
max =  1672

Ahora, si queremos muestrear los vértices de la matriz 3×3 tenemos que usar un buffer 2×2 = 4. Por tanto, hemos de sustituir las líneas 13 y 14 del código anterior por las siguientes:

.
.
.
scanline = band.ReadRaster(5, 7, 3, 3,2, 2, band.DataType)
values = struct.unpack(fmttypes[BandType] * 4, scanline)
.
.
.

La ejecución nuevamente del código resulta ahora en:

Python Console 
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(u'C:/Users/zeito/pyqgis_scripts/scanline.py'.encode('mbcs'))
(1669, 1661, 1672, 1669)
mean =  1667.75
median =  1669.0
min =  1661
max =  1672

Observen que los valores de la “cruz central” de la matriz han sido omitidos.

Esta entrada fue publicada en GDAL, PyQGIS, QGIS, 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