lunes, 20 de febrero de 2017

Rotación de vectoriales con processing grass:v.transform en PyQGIS

El comando grass:v.transform de processing permite, mediante una transformación afin, realizar translaciones y rotaciones de archivos vectoriales; independientemente de su tipo y número de rasgos. Las translaciones pueden implementarse de una manera directa pero las rotaciones involucran un grado de complejidad adicional. Ello se debe a que las rotaciones se realizan con relación al eje (0,0) y para implementarlas con relación a otro eje (por ejemplo el centroide de un vectorial) es necesario un cambio del sistema vectorial de referencia.



La matriz a continuación:



representa la rotación de θ grados del plano en sentido antihorario. Por tanto, puede demostrarse que cualquier rotación antihoraria en grados sexagesimales viene dada por:

xshift = x1 – Cos θ. x1 + Sen θ. y1

yshift = y1 – Sen θ. x1 – Cos θ. y1

donde (x1, y1) son las coordenadas del nuevo centro de referencia. Los parámetros a usar en grass:v.transform corresponden a:


El código desarrollado para una rotación arbitraria de 60 º de la capa activa corresponde a:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import processing
import math

layer = iface.activeLayer()

extent = layer.extent()

xmin, ymin, xmax, ymax = extent.toRectF().getCoords()

#rotation axis
pt = extent.center()

angle = 60

xshift = pt[0] - math.cos(angle*math.pi/180)* pt[0] + math.sin(angle*math.pi/180)* pt[1]
yshift = pt[1] - math.sin(angle*math.pi/180)* pt[0] - math.cos(angle*math.pi/180)* pt[1]

path = processing.runalg('grass:v.transform', 
                         layer,      #input <ParameterRaster>
                         xshift,           #xshift <ParameterNumber>
                         yshift,           #yshift <ParameterNumber>
                         0,                #zshift <ParameterNumber>
                         1,                #xscale <ParameterNumber>
                         1,                #yscale <ParameterNumber>
                         1,                #zscale <ParameterNumber>
                         angle,            #zrot <ParameterNumber>
                         "%f,%f,%f,%f" % (xmin, xmax, ymin, ymax), #GRASS_REGION_PARAMETER <ParameterExtent>
                         -1,               #GRASS_SNAP_TOLERANCE_PARAMETER <ParameterNumber>
                         0,                #GRASS_MIN_AREA_PARAMETER <ParameterNumber>
                         0,                #GRASS_OUTPUT_TYPE_PARAMETER <ParameterSelection> 3-area
                         None)             #OUTPUT <OutputRaster>

rotated_layer = QgsVectorLayer(path['output'],
                               'rotated',
                               'ogr')

QgsMapLayerRegistry.instance().addMapLayer(rotated_layer)

Como el comando no sobrescribe vectoriales con el mismo nombre se optó por crear uno temporal en memoria.

La ejecución del código anterior produce los resultados esperados para diferentes tipos de vectorial. En el caso del vectorial tipo polígono se tiene:


para el vectorial de puntos:


y en el caso del vectorial de línea se obtiene:



No hay comentarios: