viernes, 1 de septiembre de 2017

Generando y rotando rectángulos a partir de puntos con atributos largo, ancho y ángulo de rotación mediante PyQGIS

Cuando se quiere generar y rotar rectángulos con base a parámetros que se encuentran en la tabla de atributos de un vectorial de puntos una buena opción es mediante un script de PyQGIS. El ancho y el largo serán usados para producir un objecto de la clase QgsRectangle que será convertido en geometría y rotado con respecto al punto (que también actúa como centroide del rectángulo).

Para ejecutar el procedimiento se consideró el shapefile de la imagen siguiente donde también se puede observar la tabla de atributos con los campos largo, ancho y ángulo de rotación ya referidos.


El código completo se encuentra a continuación:

 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
38
39
40
41
42
layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

attr = [ feat.attributes() for feat in feats ]

lengths = [ attr[i][1] for i in range(len(attr)) ]
widths = [ attr[i][2] for i in range(len(attr)) ]
rotations = [ attr[i][3] for i in range(len(attr)) ]

points = [ feat.geometry().asPoint() for feat in feats ]

max_points = [ (point[0] - lengths[i]/2, point[1] + widths[i]/2) 
               for i, point in enumerate(points) ]

min_points = [ (point[0] + lengths[i]/2, point[1] - widths[i]/2) 
               for  i, point in enumerate(points) ]

rectangles = [ QgsRectangle(QgsPoint(max_points[i][0], max_points[i][1]), QgsPoint(min_points[i][0], min_points[i][1])).asWktPolygon()
               for i in range(len(points)) ]

epsg = layer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'rotated_rectangle',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(rectangles)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    geom = QgsGeometry.fromWkt(rectangles[i])
    geom.rotate(rotations[i], points[i])
    feat.setGeometry(geom)

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

La ejecución del código anterior en la Python Console de QGIS produce, tal como se esperaba, el resultado siguiente:


No hay comentarios: