Medias por columna de un ráster usando GDAL en ambiente de PyQGIS

En un artículo anterior se consideró el cálculo de las medias por filas de un ráster, en ambiente de PyQGIS, empleando los módulos GDAL y struct de Python para acceder a los valores del arreglo matricial. El cálculo de medias por fila se hizo fácil porque esa es la orientación de barrido mediante la scanline. Esta emplea la función ‘ReadRaster’ cuya sintaxis, en el caso anterior, fue la siguiente:

scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType)

donde los argumentos corresponden, respectivamente, al origen 0 en X (es decir la columna 0), el origen Y en Y (filas 0, 1, 2, …, ya que estamos en un bucle), tamaño de bloque X, tamaño de bloque Y, buffer del tamaño de bloque X, buffer del tamaño de bloque Y, el tipo de datos.

Si queremos medias por columnas entonces la scanline debe barrer según esta orientación y el código en el bucle cambiaría a:

scanline = band.ReadRaster(x, 0, 1, band.YSize, 1, band.YSize, band.DataType)

La scanline de retorno es de tipo string y contiene XSize*4 bytes de datos binarios en punto flotante. Por medio de struct.unpack se transforman estos largos strings en floats (pero ahora mediante band.YSize en lugar de band.XSize) y se almacenan en el arreglo values.

El código completo es el siguiente:

from osgeo import gdal
import os
import struct
 
layer = iface.activeLayer()
provider = layer.dataProvider()
 
fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}
 
path= provider.dataSourceUri()
 
(raiz, filename) = os.path.split(path)
 
dataset = gdal.Open(path)
 
band = dataset.GetRasterBand(1)
   
totHeight = 0
totColumns = 0
 
print "filas = %d columnas = %d" % (band.YSize, band.XSize)
 
BandType = gdal.GetDataTypeName(band.DataType)
   
print "Tipo de datos = ", BandType
  
print "Ejecutando estadisticas de %s" % filename
print "que se encuentra en %s" % raiz
 
medias_columna = []
 
for x in range(band.XSize):
   
    scanline = band.ReadRaster(x, 0, 1, band.YSize,1, band.YSize, band.DataType)
    values = struct.unpack(fmttypes[BandType] * band.YSize, scanline)
   
    for value in values:
        totHeight += value
        totColumns += value
           
    medias_columna.append(totColumns/float(band.YSize))
    totColumns = 0
     
average = totHeight / float((band.XSize * band.YSize))
(min,max) = band.ComputeRasterMinMax(1)
   
print "Promedio = %0.5f" % average
print "valor minimo = %.3f  valor maximo = %.3f" % (min,max)
  
dataset = None

Su ejecución en la Python Console para un raster GraySingleBand cargado en la Map View de QGIS produce esta salida (incluyendo un slice de las primeras 11 medias por columnas):

Python Console 
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(u'C:/Users/zeito/pyqgis_scripts/read_gdal4.py'.encode('mbcs'))
filas = 791 columnas = 1680
Tipo de datos =  Int16
Ejecutando estadisticas de utah_demUTM2.tif
que se encuentra en C:/pyqgis_data
Promedio = 1824.71801
valor minimo = 1323.000  valor maximo = 3556.000
medias_columna[0:10]
[1685.0429835651075, 1681.7901390644754, 1679.0366624525916, 
1676.3577749683943, 1673.558786346397, 1670.6738305941847, 
1667.9936788874843, 1665.8103666245258, 1664.2718078381795, 
1663.3400758533503]

Una manera directa de comprobar que el código parece producir las medias esperadas por columnas es que el promedio de 1824.71801 es idéntico al obtenido con el código para determinar las medias por filas. Sin embargo, no es suficiente. Una verificación adicional cruzada consiste en almacenar los values en un arreglo matricial y calcular las medias, por ejemplo, con numpy.

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

Una respuesta a Medias por columna de un ráster usando GDAL en ambiente de PyQGIS

  1. Pingback: Subregiones ráster usando GDAL en ambiente de PyQGIS | El Blog de José Guerrero

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