Delimitación de áreas de drenaje, mediante un script de python en GRASS, con base en los sumideros

En un artículo anterior se intentó delimitar las áreas de drenaje de una zona combinando en un sólo ráster todos los individuales para la salida (con values secuenciales; flag n) de r.drain. Como los ráster individuales abarcan toda la zona de GRASS, pero sólo tienen valores distintos de null en las líneas de flujo, es conveniente convertir los null en ceros para facilitar el álgebra de mapas en su incorporación en un ráster de valor constante (-1) que abarca toda el área de trabajo de GRASS. El producto final refleja, además del área barrida, los diferentes sumideros hacia los cuales se dirige el flujo. No obstante la estrategia seguida, es difícil la delimitación precisa de las áreas de drenaje.

Sin embargo, al observar el ráster resultado, donde quedan claramente establecidos los sumideros, se intuye que si existe la manera de obtener las coordenadas de estos lugares, donde un grupo de celdas tienen un punto de encuentro común, está sería la manera ideal de delimitar las zonas de drenaje. Después de probar varias estrategias, la que resultó más prometedora fue aquella en la cual el ráster de líneas de flujo se convirtió en un vectorial de puntos donde la mayor categoría correspondía al sumidero. Con v.extract se convirtió el vectorial multipunto en un vectorial de un sólo elemento basando la extracción en la categoría máxima para luego filtrar sólo la coordenada x con v.to.db en modo print (flag p) e incorporarlo en un arreglo matricial de n filas x n columnas.

Para incorporar esos valores en un ráster, se redujo momentáneamente el área de trabajo de GRASS al área de barrido del ráster, se produjo uno temporal en formato ascii con r.out.ascii para copiar las 6 primeras líneas de su encabezado en uno de salida, también temporal, al cual se le incorporó al final el arreglo matricial de coordenadas x antes de importarlo finalmente a GRASS con r.in.ascii. El codigo python que refleja todo lo anterior es 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 *

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", 
							flags="n",
							input="N09W068_UTM19N_canoa", 
							output=map, 
							coordinate=p1,
							quiet="true",
							overwrite="true")

		grass.run_command("r.to.vect", 
							flags="v",
							input=map, 
							output=map, 
							feature="point",
							quiet="true",
							overwrite="true")

		categorias=grass.read_command("r.category",
										map=map)

		categorias=categorias.split("\n")

		cat_max=len(categorias)-1

		grass.run_command("v.extract", 
							input=map, 
							output="temporal", 
							type="point",
							list=cat_max,
							overwrite="true",
							quiet="true")

		coord=grass.read_command("v.to.db",
							flags="p",
							map="temporal",
							type="point",
							option="coor",
							units="me",
							quiet="true")

		coord=float(coord.split("|")[1])

		ras_lambda[i][j] = coord

		grass.run_command("g.remove",
							rast=map,
							quiet="true")		

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

		grass.run_command("g.remove",
							vect="temporal",
							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="lambda2",
					quiet="true",
					overwrite="true")

pfile.close()

system ("rm " + tempfile)

Cuando se ejecuta para producir un ráster de 50 filas por 50 columnas este es el resultado; sobrepuesto al ráster con los sumideros que se produjo en un artículo anterior:

Se puede observar que hay delimitadas 7 zonas de drenaje con 5 sumideros fuera de la zona del ráster. Las áreas amarillas tienen sumideros interiores.

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

Una respuesta a Delimitación de áreas de drenaje, mediante un script de python en GRASS, con base en los sumideros

  1. Pingback: Delimitación de áreas de drenaje, mediante un script de python en GRASS, con base en los sumideros | #Geoprocessamento em Foco | Scoop.it

Deja un comentario

Introduce tus datos o haz clic en un icono para iniciar sesión:

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