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:
Manejo de arrays y matrices en python
para hacerlo más conciso y evitar el uso de g.region.
Pingback: Script de python para generar y rotar rectángulos a partir de puntos en QGIS-GRASS | ArcPY - Python | Scoop.it