Generalización de un script de python para convertir un archivo vectorial de puntos en un vectorial tipo polígono (GRASS-QGIS)

En el artículo anterior se consideró la conversión de un archivo vectorial de puntos en un vectorial tipo polígono usando un script de python y los módulos de sus librerías para GRASS (QGIS). En éste se exploran los métodos de los módulos core y vector (librerías Python-GRASS) y el manejo de string de caracteres para generalizar el procedimiento a cualquier vectorial de puntos del dataset.

El script es:

#! /usr/bin/env python
import os
import shapely
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

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, "\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

#Impresion, por verificacion, de la matriz
print "Puntos (x,y,z)"

for i in range(n):
	for j in range(3):
		print matriz[i][j],
	print 	

#Calculo del centroide

sumax = sumay = sumaz = 0.

for i in range(n):
	sumax += float(matriz[i][0])
	sumay += float(matriz[i][1])
	sumaz += float(matriz[i][2])

averagex = sumax/n
averagey = sumay/n
averagez = sumaz/n

print "\nCentroide\n", averagex, averagey, averagez

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

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

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

pfile3=open(tempfile2, 'w')  #Abre el archivo temporal

#Imprime archivo temporal para el area tipo poligono

print >> pfile2, "VERTI:\nB ", str(n)

for i in range(n):
	print >> pfile2, str(matriz[i][0]), str(matriz[i][1]), str(matriz[i][2])

print >> pfile2, "C 1 1"

print >> pfile2, averagex, averagey, averagez

print >> pfile2, "1 1"

#Imprime archivo temporal para el centroide

print >> pfile3, "VERTI:\nP 1 1"

print >> pfile3, averagex, averagey, averagez

print >> pfile3, "1 1"

pfile2.close() #Cierra el archivo temporal
pfile3.close() #Cierra el archivo temporal

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

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

print "El vectorial tipo centroide es = ", vectorial_centroide

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

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

#remueve los temporales 
os.system("rm " + tempfile)

os.system("rm " + tempfile2)

y su ejecución en cónsola de GRASS-QGIS produce la siguiente salida:

y el despliegue de los siguientes vectoriales:

Nota: El script, con pocas modificaciones, se corrió en Windows. Sólo produjo el shapefile con las coordenadas del centroide a pesar de que el archivo temporal del tipo polígono tenía el formato adecuado. Por ello, se intentó cargarlo directamente, en cónsola de GRASS, con v.in.ascii y tampoco fue posible.

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

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