Calcular media aritmética por columnas en un ráster usando un script de python

Hace algún tiempo desarrollé un script de bash para calcular la media aritmética por columnas en un ráster en ambiente de GRASS-QGIS. El script, desde el SIG y producto más bien de un reto intelectual, tardaba casi 16 horas procesando los datos lo cual es un tiempo exageradamente largo si se compara con la ejecución expedita empleando QGIS y lenguaje R. Como hay algo que se denomina criterio, es obvio que nadie en su sano juicio va a preferir usar un procedimiento de más de 15 horas, así se ejecute desde el SIG, en perjuicio de hacerlo externamente en pocos minutos. Sin embargo, existe la posibilidad de hacerlo dentro del entorno de GRASS-QGIS (aunque este no es el caso) en pocos segundos si se emplea un simple script de python.

El ráster *.tif de los dos artículos anteriores se convirtió en *.asc en QGIS y se abrió con gedit para verificar que el encabezado de la metadata del archivo tenía 5 líneas (razón del por qué un skip_header=5), 3630 filas y 3603 columnas. Se empleó en el siguiente código python:

#! /usr/bin/env python

#Librerias
import numpy as np
from StringIO import StringIO 
from scipy.stats import *

pfile1=open('tiznados_canoa.asc','r')

print "Leyendo datos..."
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..."

pfile2=open('tiznados_medias_col','w')

for i in range (3603):
	data2 = data[range(3630),i]
	media=data2.mean()
	pfile2.write("%.8f" % media)
	pfile2.write("\n")

pfile2.close()

Su ejecución en cónsola de bash, en poco más de 20 segundos, produjo la siguiente salida:

zeito@debian:~$ time python py_data3.py
Leyendo datos...
Ejecutando estadisticas...

real	0m20.058s
user	0m18.721s
sys	0m0.680s

donde el tiznados_medias_col que albergaba las medias por columnas del ráster era idéntico al obtenido con los procedimientos de los dos artículos anteriores.

El código de python no sólo es más conciso y legible que el de bash sino extremadamente eficiente en cuanto al tiempo de ejecución.

Esta entrada fue publicada en Código Python, SIG, Teledetección. Guarda el enlace permanente.

Una respuesta a Calcular media aritmética por columnas en un ráster usando un script de python

  1. Pingback: Determinación de los valores promedios de un ráster usando el módulo GDAL de Python |

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