Determinación de los valores promedios de un ráster usando el módulo GDAL de Python

En un artículo precedente se calculó la media aritmética por columnas en un ráster, en formato ascii, usando un script de python. En el presente, se va a usar el código de la página 51 de Python Geospatial Development, ligeramente modificado para adaptarlo a mis condiciones locales, para determinar el promedio total del mismo ráster. Este código usa, primero, el módulo gdal de Python para obtener y leer la banda simple del ráster a través de una línea de barrido y luego el módulo struct para leer los valores individuales fuera de la línea referida. La ventaja, con relación al primer código, es que permite usar directamente el ráster en su formato *.tif original; sin necesidad de convertirlo previamente a ascii. Además, tiene métodos que permiten extraer el número de filas y columnas del ráster y su tiempo de ejecución es considerablemente menor.

Con fines de comparación se va a modificar el código que no usa gdal ni struct para que produzca sólo la suma total de píxeles del ráster (un DEM) y su promedio. El código es el siguiente:

#! /usr/bin/env python
# -*- coding: utf-8

#Librerias
import numpy as np
from StringIO import StringIO 
from os import system

system("clear")

nameraster = 'tiznados_canoa.asc'

pfile1=open(nameraster,'r')

print "Leyendo datos de %s ..." % nameraster
data=pfile1.read()

pfile1.close()

#Se sobre entiende que los delimitadores son espacios
#El encabezado del archivo a desechar son 5 lineas 

data=np.genfromtxt(StringIO(data), skip_header=5)  
                                                  
print "Ejecutando estadisticas..."

sum = 0

for i in range (3630):
	for j in range (3603):
		sum += data[i][j]

print "Suma de píxeles = %d  Promedio = %0.5f" % (sum, sum/(3630*3603))

cuya ejecución, en cónsola de bash, produce al cabo de unos 30 segundos lo siguiente:

Leyendo datos de tiznados_canoa.asc ...
Ejecutando estadisticas...
Suma de píxeles = 3353546486  Promedio = 256.40911

real	0m30.104s
user	0m29.498s
sys	0m0.564s

El código de la página 51 de Python Geospatial Development, ligeramente modificado, es el siguiente:

#!/usr/bin/env python
# -*- coding: utf-8

from osgeo import gdal
from os import system
import struct

system("clear")

nameraster = "tiznados.tif"

dataset = gdal.Open(nameraster)
band = dataset.GetRasterBand(1)

fmt = "<" + ("h" * band.XSize)

totHeight = 0

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(fmt, scanline)

    for value in values:
        totHeight += value
		
average = totHeight / float((band.XSize * band.YSize))

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

print "Suma de píxeles = %d  Promedio = %0.5f" % (totHeight, average)

y su ejecución en cónsola de bash produce la siguiente salida:

Ejecutando estadisticas de tiznados.tif
columnas = 3603 filas = 3630
Suma de píxeles = 3353546486  Promedio = 256.40911

real	0m2.793s
user	0m2.768s
sys	0m0.024s

donde se observa que los resultados son idénticos a los del primer código pero se ejecuta en un tiempo significativamente menor y de casi 3 segundos.

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

2 respuestas a Determinación de los valores promedios de un ráster usando el módulo GDAL de Python

  1. Pingback: Cómo leer los valores de un ráster usando el módulo GDAL de Python |

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