Script de Python-GRASS-QGIS para obtener las coordenadas de los extremos de una línea y crear una nueva capa de puntos a partir de ellas

El amigo gene de GeoTux me señaló las bondades de la integración Python-GRASS recomendándome la lectura del artículo siguiente:

GRASS GIS v. 6.4.2: la nouvelle console Python (Shell Python), exemple d’utilisation simple ou avec le module Shapely

Aunque el mismo está en francés el lenguaje es sencillo por lo que se entiende fácilmente (lo difícil es lo técnico :)). Refiere la posibilidad de obtener las coordenadas de los extremos de una línea y crear una nueva capa de puntos a partir de ellas usando la nueva cónsola de Python que acompaña a la versión 6.4.2 de GRASS GIS. Como uso la versión estable de Debian mi versión de GRASS es inferior a ésta pero tampoco estoy interesado en salir de mi entorno de GRASS-QGIS. Por ello, adapté el procedimiento del artículo al ambiente de este último; ya que desde la cónsola de GRASS-QGIS puedo ejecutar el intérprete de Python (v. 2.6.6-3) o scripts códificados en este lenguage.

Para resumir el problema, lo que se pretende aquí es extraer las coordenadas de principio y fin de una línea mediante el filtrado de la salida de v.to.db con la función miembro split, crear un archivo temporal con el formato que soporta v.in.ascii en modo estándar de GRASS y, finalmente, añadirlos a una capa de puntos creada y editada mediante v.edit. A diferencia del artículo en francés, donde se trabaja en modo interactivo a través de cónsola de Python, aquí se propone el siguiente script que se ejecuta en la cónsola de GRASS-QGIS (bash) mediante python py_prueba.py (py_prueba.py es el nombre con el que decidí grabar el script y por el prefijo py_ los scripts de python estarán más o menos ordenados en el mismo sitio en mi home de usuario).

#! /usr/bin/python
import os
import grass.script as grass

os.system("clear")

#coordenadas del primer punto
startline =grass.read_command("v.to.db", flags="p", map="testgrass", type="line", option="start", units="meters" , quiet=True)

print "Esta es la salida sin filtrar de v.to.db = startline\n"

print startline

startx=float(startline.split("|")[1])
starty=float(startline.split("|")[2])
startz=float(startline.split("|")[3])

print "Despues de ejecutarse startline.split\n"

print startx, starty, startz

#coordenadas del segundo punto
endline =grass.read_command("v.to.db", flags="p", map="testgrass", type="line", option="end", units="meters" , quiet=True)

print "\nEsta es la salida sin filtrar de v.to.db = endline\n" 

print endline

endx=float(endline.split("|")[1])
endy=float(endline.split("|")[2])
endz=float(endline.split("|")[3])

print "Despues de ejecutarse endline.split\n"

print endx, endy, endz

vectorial = raw_input("\nNombre del vectorial a crear = ? ")

grass.read_command("v.edit", tool="create", map=vectorial, quiet="True")

#Creacion de un archivo temporal con las coordenadas de los puntos
tempfile = grass.tempfile()

pfile=open(tempfile, 'w')

print >> pfile, 'P', '1', '1'
print >> pfile, startx, starty, startz
print >> pfile, '1', '1'
print >> pfile, 'P', '1', '1'
print >> pfile, endx, endy, endz
print >> pfile, '1', '2'

pfile.close()

grass.run_command("v.edit", input=tempfile, flags="n", tool="add", map=vectorial, snap="node", thresh="1", quiet="True")

Para usar la misma nomenclatura que en el artículo francés, lo apliqué a un vectorial tipo línea de nombre testgrass para producir la capa vectorial de puntos de nombre testpt (aunque aquí el script me permitía cualquier otro nombre). El resultado se encuentra a continuación:

La imagen que se observa de la cónsola de GRASS-GIS es para visualizar los efectos de split sobre las salidas de v.to.db. Esos print pueden eliminarse sin problemas del script.

Anuncios