Script de python para generar y rotar rectángulos a partir de puntos en QGIS-GRASS

El siguiente código python es una traducción del correspondiente a bash basado en el mismo objetivo: generar y rotar rectángulos a partir de puntos en QGIS-GRASS. La motivación fue producto de la curiosidad en conocer la manera de traducir la instrucción al lenguaje de grass.script correspondiente a db.execute. Ello es factible, a modo de ejemplo, con algo similar a esto:

grass.write_command("db.execute", stdin = 'UPDATE rectangulo_area SET x=coor_x WHERE cat = 1')

Sin embargo, fue imposible lograr que x tomara el valor de coor_x; a menos que lo colocara explícitamente en la instrucción. Finalmente opté por v.db.update. Este es el script:

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

import os
from math import *
import grass.script as grass
import numpy as np
from StringIO import StringIO

os.system("clear")

grass.use_temp_region() #Al finalizar el script recupera la region original de grass

#El comando v.to.db se ejecuta para determinar coordenadas en modo impresion
#No requiere campos x,y en tabla atributiva

n = int(raw_input("Número de puntos pi a procesar = ? "))

print "Espere..."

for i in range(n):

	name_map = "p"+str(i+1)

	puntos=grass.read_command("v.to.db", 
								flags="p", 
								map=name_map, 
								type="point", 
								option="coor", 
								units="me", 
								quiet="True")

	puntos=np.genfromtxt(StringIO(puntos), delimiter="|")

	coor_x=puntos[1]
	coor_y=puntos[2]

	#El comando v.db.select se ejecuta para determinar el valor del radio y ang_rot

	temp=grass.read_command("v.db.select", 
								map=name_map, 
								columns="d,ang_rot")

	temp=np.genfromtxt(StringIO(temp), delimiter="|", skip_header=1)

	d=temp[0]
	ang_rot=temp[1]

	grass.run_command("v.transform", 
						input=name_map, 
						output="V1", 
						xshift="-2*d", 
						yshift="4*d", 
						quiet="True", 
						overwrite="True")

	grass.run_command("v.transform", 
						input=name_map, 
						output="V2", 
						xshift="2*d", 
						yshift="4*d", 
						quiet="True", 
						overwrite="True")

	grass.run_command("v.transform", 
						input=name_map, 
						output="V3", 
						xshift="2*d", 
						yshift="-4*d", 
						quiet="True", 
						overwrite="True")

	grass.run_command("v.transform", 
						input=name_map, 
						output="V4", 
						xshift="-2*d", 
						yshift="-4*d", 
						quiet="True", 
						overwrite="True")

	grass.run_command("v.patch", 
						input="V1,V2,V3,V4", 
						output="rectangulo_puntos", 
						quiet="True", 
						overwrite="True")

	grass.run_command("g.remove", 
						vect="V1,V2,V3,V4", 
						quiet="True")

	grass.run_command("g.region", 
						vect="rectangulo_puntos")

	grass.run_command("v.in.region", 
						type="area", 
						output="rectangulo_area", 
						quiet="True", 
						overwrite="True")

	grass.run_command("g.remove", 
						vect="rectangulo_puntos", 
						quiet="True")

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

	#Agrega las columnas x,y,ang_rot,d al vectorial tipo poligono

	grass.run_command("v.db.addcol", 
						map="rectangulo_area", 
						columns="x double precision,y double precision,ang_rot double precision,d double precision", 
						quiet="True")

	grass.run_command("v.db.update", 
						map="rectangulo_area", 
						col="x", 
						qcol=coor_x)

	grass.run_command("v.db.update", 
						map="rectangulo_area", 
						col="y", 
						qcol=coor_y)

	grass.run_command("v.db.update", 
						map="rectangulo_area", 
						col="d", qcol=d)

	grass.run_command("v.db.update", 
						map="rectangulo_area", 
						col="ang_rot",
						qcol=ang_rot)

	#Rotacion del vectorial tipo poligono

	pi=4*atan(1)

	ang_rot = 360 - ang_rot

	a= coor_x - cos(ang_rot*pi/180)*coor_x + sin(ang_rot*pi/180)*coor_y

	b= coor_y - sin(ang_rot*pi/180)*coor_x - cos(ang_rot*pi/180)*coor_y

	grass.run_command("v.transform",
						input="rectangulo_area",
						output="rectangulo_area_rot"+str(i+1),
						xshift=a, 
						yshift=b,
						zrot=ang_rot,
						quiet="True",
						overwrite="True")

	grass.run_command("g.remove", 
						vect="rectangulo_area", 
						quiet="True")

print "Listo!"

cuya ejecución para 5 puntos (p1, p2, p3, p4 y p5) produjo el siguiente resultado:

En la imagen superior los puntos aparecen señalados en sentido contrario a las agujas del reloj donde también se señala el ángulo de rotación correspondiente. Los objetivos de este script se puede abordar también con los recursos empleados en:

Conversion de un archivo vectorial de puntos en un vectorial tipo polígono usando un script de python (GRASS-QGIS)

Manejo de arrays y matrices en python

para hacerlo más conciso y evitar el uso de g.region.

Anuncios