miércoles, 24 de mayo de 2017

Dividiendo círculos en 12 secciones mediante PyQGIS

Para dividir un círculo en 12 doce secciones se puede recurrir a la obtención de un buffer con 12 lados, de un vectorial de punto, para luego producir cada una de las divisiones empleando las coordenadas de los puntos exteriores y las de su centroide (este último común a todas ellas). Como la visualización de la capa resultante no es enteramente circular, paralelamente mediante el método buffer de QgsGeometry, se produce un buffer algo menor para que al ser intersectado con cada rasgo 1/12 se obtenga el efecto deseado.


El código completo se incluye 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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import psycopg2
import numpy as np

bufferLength = 600
polygonSides = 12

registry = QgsMapLayerRegistry.instance()

layer = registry.mapLayersByName('point')

feat_points = [ feat for feat in layer[0].getFeatures() ] 
points = [ feat.geometry().asPoint() for feat in layer[0].getFeatures() ]

epsg = layer[0].crs().postgisSrid()

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

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

prov = mem_layer.dataProvider()

for i, point in enumerate(points):
    outFeat = QgsFeature()

    outFeat.setGeometry(QgsGeometry.fromPolygon([[ QgsPoint(point[0] + np.sin(angle)*bufferLength, point[1] + np.cos(angle)*bufferLength)
                        for angle in np.linspace(0, 2*np.pi, polygonSides, endpoint = False) ]]))

    outFeat.setAttributes([i])
    prov.addFeatures([outFeat])

feats_mem = [ feat for feat in mem_layer.getFeatures() ]

mem_layer2 = QgsVectorLayer(uri,
                           'sections',
                           'memory')

prov = mem_layer2.dataProvider()

k = 0

for feat in feats_mem:
    geom = feat.geometry().asPolygon()[0]
    n = len(geom)
    new_pol = []
    for i in range(n-1):
        new_pol.append([[ points[k], geom[i], geom[i+1]]])

    feat = QgsFeature()

    buffer_geom = feat_points[k].geometry().buffer(500, -1)

    for i, element in enumerate(new_pol):
        feat = QgsFeature()
        geom = QgsGeometry.fromPolygon(element)
        new_geom = geom.intersection(buffer_geom)
        feat.setGeometry(new_geom)
        feat.setAttributes([i])
        prov.addFeatures([feat])

    QgsMapLayerRegistry.instance().addMapLayer(mem_layer2)

    k += 1

Después de ejecutado en la Python Console de QGIS, para una capa vectorial de puntos con sólo tres rasgos, se obtiene lo siguiente:


No hay comentarios: