Uso del módulo Shapely (python) para determinar geometrías: comparativa frente a grass.script

A raíz de la lectura del artículo de Martin Laloux me decidí a programar en python usando las bondades del módulo grass.script. Sin embargo, en ese artículo también se refiere un uso simple del módulo Shapely para crear objetos tipo punto y línea, la obtención de algunos de sus parámetros de geometría (longitud de línea) y su representación bidimensional usando matplotlib. La ventaja que tiene el uso del módulo Shapely es que se pueden determinar geometrías asignadas a objetos tipo polígono, como área y perímetro, aprovechando la lista correspondiente a la individualización de coordenadas producto de la salida de v.to.db sin necesidad de crear esos objetos espaciales en el espacio de trabajo de GRASS.

Para instalar Shapely y matplotlib en Debian:

apt-get install python-shapely
apt-get install python-matplotlib

Para Windows XP bajé Shapely (para python 2.7.x) de aquí:

http://pypi.python.org/pypi/Shapely#downloads

y después de su instalación se cargó sin problemas en el intérprete.

El script usado para probar Shapely en Debian fue el siguiente:

#! /usr/bin/env python
import os
from shapely.geometry import Polygon
from shapely.geometry import Point
from grass.script import core as grass
from grass.script import vector as grass

os.system("clear")

lista = grass.list_strings('vect') #Lista de vectoriales en el mapaset actual
long_lista=len(lista) #longitud de la lista

j=0

for i in range(long_lista):
	info = grass.vector_info_topo(lista[i])
	if (info['lines'] == 0 and info['boundaries'] == 0):
		lista[j]=lista[i]
		j += 1

long_lista = j #nueva longitud de la lista (solo vectoriales de punto)

ver=lista[0] #Primer elemento de la lista
long_first=len(ver) #longitud primer elemento de la lista

cont = 0 #cont es la posicion de @ en el primer elemento de la lista

while(ver[cont] != '@'):
	cont += 1

subcadena_ver=ver[cont:len(ver)] #Representa el mapaset actual
long_mapaset=len(subcadena_ver)  #Longitud de la cadena del mapaset actual

print "Lista de los " + str(long_lista) + " vectoriales de puntos en el mapset actual ("+ subcadena_ver + ")\n"

h = 4

for i in range(long_lista):
	l_temp=len(lista[i])
	lista_temp=lista[i]
	lista[i]=lista_temp[0:(l_temp-long_mapaset)] #elimina de las cadenas el @mapaset	
	print lista[i],
	if (i == h):
		print
		h += 5 

print "\n"

vectorial = raw_input ("Nombre del vectorial a procesar = ? ")

info = grass.vector_info_topo(vectorial) #Diccionario con la info del vectorial 
                                         #Numero de puntos es 'points'

#numero de puntos

n = info['points']

print "\nEl numero de puntos del vectorial " + vectorial + " es =", n

#Obtencion de las coordenadas x,y,z para todos los puntos
puntos = grass.read_command("v.to.db", flags="p", map=vectorial, type="point", option="coor", units="meters", quiet="True")

#Definicion de la matriz
matriz = []

for i in range(n):
	matriz.append([])
	for j in range(3):
		matriz[i].append(None)

#Algoritmo para la individualizacion de las coordenadas
z=1
h=0
for i in range(n):
	h += 3
	for j in range(3):
		matriz[i][j]=(puntos.split("|")[z])
		if (j==2):
			temp=(puntos.split("|")[h])
			matriz[i][j]=float(temp.split("\n")[0])
		z += 1


vectorial_area=vectorial + "_area"
vectorial_centroide =vectorial + "_centroide"

print "\nEl vectorial tipo poligono es = ", vectorial_area

for i in range(n):
	for j in range(3):
		matriz[i][j] = float(matriz[i][j])	

area_ = grass.read_command("v.to.db", flags="p", map=vectorial_area, type="boundary", option="area", units="meters", quiet="True")
perimetro_ = grass.read_command("v.to.db", flags="p", map=vectorial_area, type="boundary", option="length", units="meters", quiet="True")

print "\nUsando libreria de GRASS (funciona porque ya estan los poligonos creados)"

area = (area_.split("|")[1])
area = (area.split("\n")[0])

perimetro = (perimetro_.split("|")[1])
perimetro = (perimetro.split("\n")[0])

print "\nEl area de", vectorial_area, "es =", area, "m2"

print "\nEl perimetro de", vectorial_area, "es =", perimetro, "m"

poligono = Polygon(matriz) #Se crea el poligono con funcion de shapely.geometry

print "\nUsando Shapely"

print "\nEl area de", vectorial_area, "es =", poligono.area, "m2"
print "\nEl perimetro de", vectorial_area, "es =",  poligono.length, "m\n"

punto=range(n)
distancia=range(n)

#Se crean los puntos con funcion de shapely.geometry

for i in range(n):
	punto[i] = Point(matriz[i])

for i in range(n-1):
	distancia[i]=punto[0].distance(punto[i+1])

h=2

for i in range(n-1):
	print "distancia["+ str(i+1) + "]=", distancia[i],
	if i == h:
		print
		h +=3	

En el script se observa que en Shapely, para crear el objeto tipo polígono se le pasa como parámetro la lista completa matriz (con conversión previa a float) y para crear los objetos tipo punto las sub listas matriz[i]. Las geometrías (área, perímetro, distancia entre puntos) se determinan con instrucciones tan sencillas como:

poligono.area

poligono.length

for i in range(n-1):
	distancia[i]=punto[0].distance(punto[i+1])

sin necesidad de incorporarlos en el ambiente de trabajo de GRASS. Como el vectorial tipo polígono fue creado como objeto de GRASS con grass.script en el artículo pasado, se puede determinar su geometría mediante esta vía para compararlos con los provenientes del uso de Shapely. Los resultados de la ejecución del script se observan en la siguiente imagen:

Son esencialmente iguales los valores de área y perímetro determinados por ambas vías. Por otra parte, en la siguiente imagen se encuentran de manera explícita las distancias entre puntos referidas en el script (la última es cero porque uno de los puntos está repetido).

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

Una respuesta a Uso del módulo Shapely (python) para determinar geometrías: comparativa frente a grass.script

  1. perfecto, también presente el uso de Shapely en http://www.portailsig.org/content/python-le-module-shapely-geometries-predicats-spatiaux-analyse-spatiale-matrices-de-clementi (en francés). El autor de Shapely es Sean Gillies y su blog http://sgillies.net/blog/. Tiene una nueva librería muy interesante que se llama Fiona http://pypi.python.org/pypi/Fiona/0.8 (no funciona en Windows)
    El punto mas importante es que todo se puede hacer en el Shell Python mezclando QGIS + GRASS + Shapely y otros sin lanzar los programas.

    import sys, os, numpy, argparse
    import grass.script as grass
    import qgis .core
    import shapely
    import fiona
    import osgeo.gdal
    import osge.ogr

    y tienes el script Python total para manejar lo que quieres

    Si quieres ejemplos de Python con GRASS, también presente cosas en:
    http://osgeo-org.1560.n6.nabble.com/access-vertices-coordinates-of-a-line-from-Python-td4983417.html (en ingles)
    http://osgeo-org.1560.n6.nabble.com/Automatic-3D-geological-boreholes-representation-automate-v-extrude-from-a-table-my-solution-in-Pythn-td4978801.html (en ingles)

    y sobre Python con QGIS
    http://www.portailsig.org/content/python-creer-automatiquement-des-fichiers-styles-qml-de-qgis-partir-de-fichiers-ou-de-shapef ((en francés)
    http://www.portailsig.org/content/qgis-lancer-des-scripts-python-ou-des-commandes-shell-depuis-la-console-python-ou-avec-scrip (en francés)

    y en
    http://plugins.qgis.org/snippets/ (en ingles)

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