Script de python para determinar el índice de alargamiento de una Cuenca Hidrográfica: GRASS-QGIS

Hace algún tiempo publiqué una saga de 3 artículos que pretendían, en conjunto, determinar el índice de alargamiento de una Cuenca Hidrográfica en QGIS. Estos fueron los siguientes:

Pendiente y punto medio, por tramos, en un shapefile tipo línea (QGIS-GRASS)

Familia de rectas ortogonales que pasan por el punto medio de cada segmento individual de una polilínea

Puntos de intersección de familia de rectas ortogonales con shapefile de cuenca hidrográfica y estimación de anchura máxima perpendicular al río

A medida que fui familiarizándome con GRASS y Python, me surgió el interés de combinar los 3 procedimientos en uno sólo mediante un script de python usando las librerías de grass.script. El objetivo es determinar el índice de alargamiento de una Cuenca Hidrográfica que se define como la relación entre el largo del río y la anchura máxima perpendicular al río. Producto de la digitalización, la recta representada por el río está constituida por una serie de segmentos lineales cuya pendiente o derivada es válida para ambos extremos (o cualquier punto interior al segmento), por lo que la distancia perpendicular al río que abarque la cuenca podría evaluarse para cualquiera de ellos. Sin embargo, en aras de la representatividad, se asume el cálculo a partir de la recta que pasa por el punto medio del segmento. Para ello, se divide el río en segmentos lineales de 2 vértices, mediante v.split, el cual los genera con una tabla atributiva de una sola categoría. Sin embargo, se requiere un vectorial con múltiples categorías y esto es solventado con la modificación de las categorías a partir de una copia temporal del vectorial en formato ascii con v.out.ascii.

El archivo ascii de partes sencillas tiene un encabezado de texto de 10 líneas y a partir de allí las modificaciones de categoría sólo son necesarias a partir de la línea 17 (categoría 2) cada 4 de ellas. La individualización de las líneas de texto se logra con split y el delimitador “\n”. Después de modificadas las líneas, se copian a un archivo temporal el cual se recompone a un vectorial de partes sencillas con categorías diferentes mediante v.in.ascii con eliminación del temporal.

Al ahora vectorial multi segmentos y multi categorías se le incorporan elementos de geometría en su tabla atributiva para determinar la familia de rectas perpendiculares que pasan por el punto medio de cada segmento y se extienden más allá del shapefile que delimita la cuenca. Esta familia de rectas se intersecta con el shapefile que delimita la cuenca usando el módulo v.overlay, se hace una extracción por categoría de cada recta, mediante v.extract, se determinan sus longitudes individuales con v.to.db y se incorporan a un arreglo donde es fácil encontrar su valor máximo.

En este punto, encontrada la distancia máxima perpendicular al río (anchura_max_perp), se determina la distancia total de este último (longitud_rio) con v.to.db para ser usada en el cálculo del índice de alargamiento como sigue:

indice de alargamiento=longitud_rio/anchura_max_perp

El script propuesto es el siguiente:

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

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

system("clear")

grass.run_command("v.split",
					input="rio",
					output="rio_split",
					vertices="2",
					quiet="True",
					overwrite="True")

puntos=grass.read_command("v.out.ascii", 
							input="rio_split", 
							format="standard", 
							quiet="True")

puntos=puntos.split("\n")

n=len(puntos)#Determina numero de lineas del arreglo puntos

#inicia contadores para modificar lineas del archivo ascii
h=0
j=2

#El primer valor a modificar esta en posicion 17 cada 4 valores
for i in range(0,len(puntos)):
	if i==17+h:
		puntos[i]="1 "+str(j) #modifica la categoria del vectorial
		h+=4
		j+=1

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

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

#Imprime archivo temporal para el vectorial tipo linea

for i in range(0,len(puntos)):
	print >> pfile, str(puntos[i])

pfile.close()

grass.run_command("v.in.ascii", 
					format="standard", 
					input=tempfile, 
					out="verizona5", 
					quiet="True", 
					overwrite="True")

system ("rm " + tempfile)

grass.run_command("v.db.addtable", 
					map="verizona5", 
					quiet="True")

grass.run_command("v.db.addcol", 
					map="verizona5", 
					columns="longitud double precision, startx double precision, starty double precision, endx double precision,endy double precision, mean_x double precision, mean_y double precision, slope double precision, slope_ort double precision, b_ort double precision")

grass.run_command("v.to.db", 
					map="verizona5", 
					type="line",
					option="length",					
					col="longitud", 
					units="me",
					quiet="True")

grass.run_command("v.to.db", 
					map="verizona5", 
					type="line", 
					option="start", 
					col="startx,starty",
					quiet="True")

grass.run_command("v.to.db", 
					map="verizona5", 
					type="line", 
					option="end", 
					col="endx,endy",
					quiet="True")

grass.run_command("v.db.update",
					map="verizona5", 
					layer="1", 
					column="mean_x", 
					qcolumn="(startx+endx)/2",
					quiet="True")

grass.run_command("v.db.update",
					map="verizona5", 
					layer="1", 
					column="mean_y", 
					qcolumn="(starty+endy)/2",
					quiet="True")

grass.run_command("v.db.update", 
					map="verizona5", layer="1", 
					column="slope", 
					qcolumn="(endy-starty)/(endx-startx)",
					quiet="True")

grass.run_command("v.db.update", 
					map="verizona5", 
					layer="1", 
					column="slope_ort", 
					qcolumn="-1/slope",
					quiet="True")

grass.run_command("v.db.update", 
					map="verizona5", 
					layer="1", 
					column="b_ort", 
					qcolumn="mean_y-slope_ort*mean_x",
					quiet="True")

grass.run_command("v.db.addcol", 
					map="verizona5",
					layer=1, 
					columns="Xn double precision,Yn double precision,Xs double precision,Ys double precision")

#Seleccion de los valores de la pendiente ortogonal
slope_ort=grass.read_command("v.db.select", flags="c", map="verizona5", col="slope_ort")

slope_ort=(slope_ort.split("\n"))

for i in range(len(slope_ort)-1):
	slope_ort[i]=float(slope_ort[i])

slope_ort= slope_ort[0:(len(slope_ort)-1)]

#Seleccion de los valores del intercepto ortogonal
b_ort=grass.read_command("v.db.select", flags="c", map="verizona5", col="b_ort")

b_ort=(b_ort.split("\n"))

for i in range(len(b_ort)-1):
	b_ort[i]=float(b_ort[i])

b_ort= b_ort[0:(len(b_ort)-1)]


#Seleccion de los valores de mean_x[i]

mean_x=grass.read_command("v.db.select", flags="c", map="verizona5", col="mean_x")

mean_x=(mean_x.split("\n"))

for i in range(len(mean_x)-1):
	mean_x[i]=float(mean_x[i])

mean_x= mean_x[0:(len(mean_x)-1)]

#Calculo de Xn[i]

d=135000 #Ancho de la cuenca

Xn = range(len(mean_x))

for i in range(len(mean_x)):
	Xn[i]=d/sqrt(slope_ort[i]*slope_ort[i]+1) + mean_x[i]

#Calculo de Yn[i]

Yn = range(len(mean_x))

for i in range(len(mean_x)):
	Yn[i] = slope_ort[i]*Xn[i] + b_ort[i]

#Calculo de Xs[i]

Xs = range(len(mean_x))

for i in range(len(mean_x)):
	Xs[i]=-d/sqrt(slope_ort[i]*slope_ort[i]+1) + mean_x[i]

#Calculo de Ys[i]

Ys = range(len(mean_x))

for i in range(len(mean_x)):
	Ys[i] = slope_ort[i]*Xs[i] + b_ort[i]

for i in range(len(mean_x)):
	input1="UPDATE verizona5 SET Xn=" + str(Xn[i]) + " WHERE cat = " + str(i+1)
	grass.write_command("db.execute", stdin = input1)

	input2="UPDATE verizona5 SET Yn=" + str(Yn[i]) + " WHERE cat = " + str(i+1)
	grass.write_command("db.execute", stdin = input2)

	input3="UPDATE verizona5 SET Xs=" + str(Xs[i]) + " WHERE cat = " + str(i+1)
	grass.write_command("db.execute", stdin = input3)

	input4="UPDATE verizona5 SET Ys=" + str(Ys[i]) + " WHERE cat = " + str(i+1)
	grass.write_command("db.execute", stdin = input4)

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

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

print >> pfile, "VERTI:"

for i in range(len(mean_x)):
	print >> pfile, "L 2 1"
	print >> pfile, str(Xn[i]), str(Yn[i])
	print >> pfile, str(Xs[i]), str(Ys[i])
	print >> pfile, "1 " + str(i+1) 

pfile.close()

grass.run_command("v.in.ascii",
					input=tempfile, 
					output="familia_lineas", 
					format="standard",
					quiet="True",
					overwrite="True")

system ("rm " + tempfile)
 
grass.run_command("v.overlay",
					ainput="familia_lineas",
					atype="line", 
					binput="cuenca", 
					btype="area", 
					output="lineas_intersec", 
					operator="and", 
					olayer="0,1,0",
					quiet="True",
					overwrite="True")

#Preparando el loop para la extraccion de las lineas en familia_lineas

largo=range(len(mean_x))

h=0

for i in range(len(mean_x)):

	h +=1

	grass.run_command("v.extract", 
						input="lineas_intersec", 
						out="temporal", 
						list=h, 
						quiet="True", 
						overwrite="True")

	salida="linea" + str(i+1)

	grass.run_command("v.build.polylines", 
					input="temporal", 
					output=salida, 
					cats="first", 
					quiet="True", 
					overwrite="True")
	
	longitud=grass.read_command("v.to.db", 
								flags="p", 
								map=salida, 
								option="length", 
								quiet="true") #quiet="true" evita encabezado

	longitud=longitud.split("|")

	longitud=longitud[1].split("\n")

	longitud=float(longitud[0])

	largo[i] = longitud

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

anchura_max_perp = max(largo)

longitud_rio=grass.read_command("v.to.db", 
							flags="p", 
							map="rio", 
							option="length", 
							quiet="true") #quiet="true" evita encabezado

longitud_rio=longitud_rio.split("|")

longitud_rio=longitud_rio[1].split("\n")

longitud_rio=float(longitud_rio[0])

indice_alargamiento=longitud_rio/anchura_max_perp

print "\nIndice de alargamiento = %.2f" % indice_alargamiento

La ejecución del script produce una serie de advertencias no fatales producto de que en tiempo de ejecución, desconociendo el motivo, se buscan las descripciones en idioma español a pesar de que la sintaxis de los comandos involucrados (v.overlay y v.build.polylines) es la correcta. El resultado, en poco más de 3 segundos, es el esperado; tal como se observa en la siguiente imagen:

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

Una respuesta a Script de python para determinar el índice de alargamiento de una Cuenca Hidrográfica: GRASS-QGIS

  1. Pingback: Script de python para determinar el índice de alargamiento de una Cuenca Hidrográfica: GRASS-QGIS | Geoprocessing - Geoprocessamento | 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