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

En un post anterior se determinaron los valores promedios de un ráster usando el módulo GDAL de Python. El método de lectura del array fue el más eficiente y utiliza una “línea de barrido” (scanline) de tipo string que contiene el número de columnas del ráster multiplicado por 4 bytes de datos binarios en punto flotante y que se convierte a valores Python utilizando el módulo struct de la librería estándar. El código usado fue adaptado de la página 51 de Python Geospatial Development y recién me doy cuenta que funcionó para mi caso particular porque el DataType de mi mapa era coincidencialmente idéntico al del considerado en el ejemplo.

El “descubrimiento” comienza con el empleo del mismo código para tratar de hacer lo mismo con un mapa *.jpg. Después de varios intentos infructuosos me decanté, en primera aproximación, por el procedimiento ineficiente que refieren aquí:

Reading Raster Data with GDAL

y produje este código:

#!/usr/bin/env python
# -*- coding: utf-8
 
from osgeo import gdal
from os import system
import struct
 
system("clear")
 
nameraster = "tiznados_canoa.jpg"
 
dataset = gdal.Open(nameraster)
band = dataset.GetRasterBand(1)

print "filas = %d columnas = %d" % (band.YSize, band.XSize)

BandType = gdal.GetDataTypeName(band.DataType)

print "Tipo de datos = ", BandType

print "Ejecutando estadisticas de %s" % nameraster

data = band.ReadAsArray(0, 0, band.XSize, band.YSize)

sum = 0

for i in range(band.YSize):
	for j in range(band.XSize):
		
		sum += float(data[i][j])
		
promedio = sum/(band.XSize * band.YSize)

print "Promedio = %0.5f" % promedio

dataset = None

que calcula el promedio en algo más de 12 segundos:

filas = 3630 columnas = 3603
Tipo de datos =  Byte
Ejecutando estadisticas de tiznados_canoa.jpg
Promedio = 173.12634

real	0m12.750s
user	0m12.665s
sys	0m0.076s

No obstante, seguí investigando la naturaleza del código y descubrí que la función miembro unpack del módulo struct depende de los fmttypes asignados a cada DataType que en GDAL son los 12 siguientes:

Unknown
Byte
UInt16
Int16
UInt32
Int32
Float32
Float64
CInt16
CInt32
CFloat32
CFloat64

Después de averiguar los fmttypes de los 7 primeros (excluyendo Unknown), creé el siguiente diccionario:

    fmttypes = {‘Byte’:’B’, ‘UInt16′:’H’, ‘Int16′:’h’, ‘UInt32′:’I’, ‘Int32′:’i’, ‘Float32′:’f’, ‘Float64′:’d’}

que incluido en el código siguiente:

#!/usr/bin/env python
# -*- coding: utf-8
 
from osgeo import gdal
from os import system
import struct
 
system("clear")
 
nameraster = "tiznados_canoa.jpg"
 
dataset = gdal.Open(nameraster)
band = dataset.GetRasterBand(1)
 
totHeight = 0

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}

print "filas = %d columnas = %d" % (band.YSize, band.XSize)

BandType = gdal.GetDataTypeName(band.DataType)
 
print "Tipo de datos = ", BandType

print "Ejecutando estadisticas de %s" % nameraster
 
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

produce el mismo resultado que el código anterior pero en menos de 3 segundos:

filas = 3630 columnas = 3603
Tipo de datos =  Byte
Ejecutando estadisticas de tiznados_canoa.jpg
Promedio = 173.12634

real	0m2.167s
user	0m2.144s
sys	0m0.020s

En principio, el código debería funcionar para estos DataType: Byte, UInt16, Int16, UInt32, Int32, Float32 y Float64. Se corroboró su funcionamiento adecuado para mapas con los siguientes DataType: Byte, Int16 y Float32.

Esta entrada fue publicada en Código Python, GDAL, SIG, Software Libre. Guarda el enlace permanente.

6 respuestas a Cómo leer los valores de un ráster usando el módulo GDAL de Python

  1. Pingback: Lista con el nombre de los archivos de una determinada extensión en un directorio mediante Python |

  2. Pingback: Cómo extraer las coordenadas x, y, z de un ráster *.tif para visualización 3D con Python |

  3. Hola José, conoces alguna manera de obtener la altura Z de un raster proporcionando la posición X,Y (Long,Lat); he intentado leer la información de GDAL pero me parece algo confusa. Requiero leer archivos *.grd y el lenguaje que manejo es VB.Net.

  4. Si mal no recuerdo, eso está resuelto en mi Blog empleando Python (con GDAL) y GRASS scripting. Busca por las etiquetas SIG, GRASS, Python. Con relación a VB.net no lo conozco.

  5. Pingback: Cómo leer los valores de un ráster usando el módulo GDAL con PyQGIS | El Blog de José Guerrero

  6. 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