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

En el artículo precedente se estableció un procedimiento para determinar la longitud del terreno (λ de la RUSLE), según enfoque CALSITE, mediante r.drain en GRASS-QGIS para celdas individuales del ráster. En éste se automatiza el procedimiento, mediante un script de python, para barrer un área del ráster DEM y producir otro, en la misma región, que contenga los valores de λ.

La lógica del programa permite 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 crea una región temporal que se alinea 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 convierte internamente a ascii (en tiempo de ejecución), se le sustituyen los valores constantes por los de la lista de valores de λ a través de un archivo temporal que luego se importa al ambiente de trabajo de GRASS.

El script propuesto (py_longitud.py) se encuentra a continuación y se usó para producir un ráster de 50 filas por 50 columnas (resolución 30.61054474 m x 30.61054474 m) a partir de la coordenada x=640096.00, y=1100717.00.

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

import grass.script as grass
import numpy as np
from os import system
from math import *

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 range(rows):
	
	print "Procesando fila " + str(i+1) 	
	
	for j in range(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.report",
									flags="hn",
									map=map,
									units="me",
									nsteps="1",
									quiet="true")

		salida=salida.split("TOTAL")

		salida=salida[1].split("|")

		salida=salida[1].split(",")

		sum_texto=""

		for k in range(len(salida)):

			sum_texto=sum_texto + salida[k]

		area_acum=float(sum_texto)

		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	

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")

grass.run_command("v.in.region",
					output="vector",
					quiet="true",
					overwrite="true")

grass.run_command("v.to.rast",
					input="vector",
					output="vector_raster",
					value="1",
					use="val",
					quiet="true",
					overwrite="true")

grass.run_command("g.remove",
					vect="vector",
					quiet="true")

salida=grass.read_command("r.out.ascii", 
					input="vector_raster")

salida=salida.split("\n")

tempfile = grass.tempfile() #Define el archivo temporal

pfile=open(tempfile, 'w')  #Abre el archivo temporal

for i in range(6):

	print >> pfile, salida[i]

for i in range(rows):

	for j in range(columns):

		print >> pfile, ras_lambda[i][j],

	print >> pfile, " "

pfile.close()

pfile=open(tempfile, 'r')  #Abre el archivo temporal

grass.run_command("r.in.ascii",
					input=tempfile,
					output="lambda",
					quiet="true",
					overwrite="true")

pfile.close()

system ("rm " + tempfile)

El resultado de ejecución del script se encuentra en la siguiente imagen (225 Has):

Las zonas en azul registran valores tan bajos como 17 m por lo que asumir valores de λ de 50 o 100 m corresponderían a una sobre estimación en dichas regiones del ráster.

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

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

  1. Jaime dijo:

    Muy bueno tu blog, y me presento, soy Jaime Talavera García, Ingeniero Forestal, ye me gustaria aprender mas sobre la RUSLE, y me pregunta es se puede trabajar con imagenes de satelite?

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

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