Cómo leer los valores de un ráster usando el módulo GDAL con PyQGIS

En un artículo anterior se exploró la necesidad de obtener las estadísticas de un ráster utilizando el método ‘bandStatistics’ del objeto QgsRasterDataProvider correspondiente. Este contempla una serie de parámetros descriptores como la media, suma de cuadrados, etc para el ráster completo. No obstante, con el módulo gdal, es posible acceder a los valores individuales del ráster y proceder a efectuar estadísticas más particulares; si es nuestro deseo.

En este artículo sólo se va a señalar como acceder al array dentro de la Python Console de QGIS empleando el módulo gdal y el cálculo sencillo del promedio de los valores de sus píxeles. Para ello, necesitamos primero acceder a la ruta (path) del ráster que tengamos cargado en la Map View de QGIS y eso se logra a través del objeto QgsRasterDataProvider ahora empleando el método ‘dataSourceUri’. Esto permite encontrar la ruta hasta el archivo correspondiente. Ahora, modificando ligeramente el script que se encuentra en otro artículo anterior, se puede acceder a los los valores del ráster usando el módulo GDAL de Python desde la propia Python Console de PyQGIS.

El código completo es:

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

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
  
for y in range(band.YSize):
  
    scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType)
    values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)
  
    for value in values:
        totHeight += value
          
average = totHeight / float((band.XSize * band.YSize))
  
print "Promedio = %0.5f" % average
 
dataset = None

el cual ejecuté desde el Editor de la Python Console habiendo cargado previamente el ráster cuyo filename es tiznados_canoa.tif. El script (gdal_read.py en mi directorio de scripts señalado en el PYTHONPATH) permite determinar el valor promedio para todos los valores del ráster.

El resultado completo de la ejecución es este:

Python Console 
Use iface to access QGIS API interface or Type help(iface) for more info
execfile(u'/home/zeito/scriptspyqgis/gdal_read.py'.encode('UTF-8'))
filas = 3630 columnas = 3603
Tipo de datos =  Int16
Ejecutando estadisticas de tiznados_canoa.tif
que se encuentra en /home/zeito
Promedio = 256.40911

del cual se anexa también una imagen:

gdal_read

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

4 respuestas a Cómo leer los valores de un ráster usando el módulo GDAL con PyQGIS

  1. Pingback: Medias por fila de un ráster usando GDAL en ambiente de PyQGIS | El Blog de José Guerrero

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

  3. Pingback: Creando un plugin de QGIS para determinar el promedio de un ráster con GDAL | El Blog de José Guerrero

  4. Pingback: Cómo producir tiles idénticos de un raśter de forma masiva con GDAL y 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