Cálculo de la longitud del terreno (λ de la RUSLE) según enfoque CALSITE mediante un script de Python con GRASS y GDAL

Hace algo más de un año establecí un procedimiento para determinar la longitud del terreno (λ de la RUSLE), según enfoque CALSITE, mediante r.drain en GRASS-QGIS. La lógica del programa permitía barrer el ráster DEM en sentido horizontal y vertical para determinar individualmente cada valor de λ y asignarlos a una lista (array de arrays). Posteriormente, se creaba una región temporal que se alineaba con el ráster DEM con el fin de obtener, primero, un vectorial de la zona barrida y luego un ráster de valor constante con la misma resolución del DEM. Este ráster se convertía internamente a ascii (en tiempo de ejecución), se le sustituían los valores constantes por los de la lista de valores de λ a través de un archivo temporal que luego se importaba al ambiente de trabajo de GRASS.

Sin embargo, con el uso de GDAL, al tener la lista de valores de λ (array de arrays) ya no es necesaria el álgebra de mapas a través de un ráster ascii de valores constantes que se sustituían por los de λ a través de un archivo creado temporalmente y que luego se importaba al ambiente de trabajo de GRASS. Con GDAL se crea un ráster *.tif (en la dirección de nuestra preferencia y fuera del ambiente de GRASS) al cual directamente se le escribe el array en una sola instrucción:

.
.
.
dst_ds.GetRasterBand(1).WriteArray( ras_lambda )
.
.
.

Posteriormente, se copia su extensión, el ancho de celda y la referencia espacial. El código fue el siguiente:

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

import grass.script as grass
import numpy as np
from os import system
from math import *
import time
from osgeo import osr, gdal

system("clear")

x=640096.00
y=1100717.00

rns=30.61054474
rew=30.61054474

n = y
w = x

rows=int(raw_input("número de filas = ? "))
columns=int(raw_input("número de columnas = ? "))

ras_lambda=np.zeros((rows,columns))

print "Espere...\n"

for i in xrange(rows):
	
	start = time.time()

	print "Procesando fila " + str(i+1)  	
	
	for j in xrange(columns):

		p1=str(x)+","+str(y)

		map="flujo"+str(i)+str(j)

		grass.run_command("r.drain", 
							input="N09W068_UTM19N_canoa", 
							output=map, 
							coordinate=p1,
							quiet="true",
							overwrite="true")

		salida = grass.read_command("r.sum",
									rast=map,
									quiet="true")

		salida = salida.split('=')

		area_acum = (float(salida[1]))*rns*rew

		pi=4*atan(1)

		ras_lambda[i][j]=sqrt(area_acum/pi)

		grass.run_command("g.remove",
							rast=map,
							quiet="true")		
		x += rew
	
	x=w
	y -= rns

	end = time.time()

	print "tiempo de proceso = %.2f minutos" % ((end - start)/60.)	

grass.use_temp_region()

s = n - (rows-1)*rns
e = w + (columns-1)*rew

grass.run_command("g.region", 
					n=n, 
					w=w, 
					s=s, 
					e=e, 
					align="N09W068_UTM19N_canoa")

lista = grass.read_command("g.region",
					flags="g")

lista = lista.split('\n')

norte = lista[0].split('=')
oeste = lista[2].split('=')

norte[1] = float(norte[1])
oeste[1] = float(oeste[1])

# Establezca archivo de salida
output_file = "/home/zeito/Desktop/raster/salida.tif"
 
# Create gtif
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, columns, rows, 1, gdal.GDT_Float64)

dst_ds.GetRasterBand(1).WriteArray( ras_lambda )

# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform( [ oeste[1], rew, 0, norte[1], 0, -rns ] )

# set the reference info 
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS( 'WGS84' )
srs.SetUTM( 19, 1 )
dst_ds.SetProjection( srs.ExportToWkt() )

resultando este mapa:

lambda

que concuerda con el obtenido anteriormente y con un tiempo de ejecución aproximado de 8 horas con 20 minutos. El objetivo es ahora optimizar el tiempo de ejecución usando otras opciones de r.drain.

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

4 respuestas a Cálculo de la longitud del terreno (λ de la RUSLE) según enfoque CALSITE mediante un script de Python con GRASS y GDAL

  1. rengifo dijo:

    Excelente artículo. Solo tengo una pregunta. Cuando usas el símbolo λ me imagino que te refieres a la longitud de la pendiente.

    • Gracias por tu comentario. En cuanto a λ es correcta tu apreciación, es decir, es la longitud de la pendiente en metros. No obstante, yo uso en el post la denominación que empleó Barrios (2000) en su artículo que la refiere como longitud del terreno; la cual está relacionada con el subfactor adimensional longitud del terreno L de la RUSLE.

      Saludos

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

  3. Pingback: Interpolación spline cubica con gdalwarp de mapas de precipitación obtenidos por interpolación TIN | 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