sábado, 25 de marzo de 2017

Corrigiendo geometrías con PyQGIS: segunda parte

En el post anterior se consideró la corrección de una geometría con problemas en el sentido de que no se podía generar el poligono correspondiente a partir de los rasgos de línea. Se determinó la existencia de varios gaps que fueron corregidos mediante v.clean de la Processing Tool Box. Sin embargo, esto no era la causa de los problemas.



Estos radicaban en el hecho de que los rasgos no estaban ordenados secuencialmente y, cuando esto último se realizó, todos permanecieron individualmente ordenados en el sentido contrario a las agujas del reloj (equivalente a señalar que estaban pegados "cabeza-cola" en el sentido equivocado, es decir, había que "girarlos" todos 180 grados).

El ordenamiento secuencial de los rasgos se resolvió con este código:

 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
import itertools

layer = iface.activeLayer()

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

selected_feats = []

selected_feats.append(feats[0])

del feats[0]

k = 0

while len(feats) != 0:
    for i, feat in enumerate(feats):
        d = selected_feats[k].geometry().distance(feat.geometry())
        if d == 0:
            selected_feats.append(feat)
            k += 1
            del feats[i]

lines = [ feat.geometry().asPolyline() for feat in selected_feats ]

epsg = layer.crs().postgisSrid()
    
uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           '4testing_feats_ordered',
                           'memory')

prov = mem_layer.dataProvider()

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

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(lines[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

que al ejecutarlo en la Python Console de QGIS produce el resultado esperado; tal como se observa en la imagen siguiente:


Para corregir el ordenamiento de los puntos en cada rasgo y hacer que todos siguieran la misma dirección en el sentido de las agujas del reloj, produciendo en el mismo procedimiento el polígono deseado, se usó este otro código:

 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
registry = QgsMapLayerRegistry.instance()

layer = registry.mapLayersByName("4testing_feats_ordered") 

lines = [ feat.geometry().asPolyline() for feat in layer[0].getFeatures() ]

sum = 0

for i in range(len(lines[0])-2):
    det =  (lines[0][i+1][0] - lines[0][i][0])*(lines[0][i+2][1] - lines[0][i][1])-(lines[0][i+1][1] - lines[0][i][1])*(lines[0][i+2][0] - lines[0][i][0])
    sum += det

if sum > 0:
    
    lines[0].reverse()

new_points = []

for point in lines[0]:
    new_points.append(point)

for i in range(1, len(lines), 1):
    if lines[i][-1] == new_points[-1]:
        lines[i].reverse()

    for point in lines[i]:
        new_points.append(point)

polygon = [[new_points]]

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

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

prov = mem_layer.dataProvider()

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

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolygon(polygon[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Después de ejecutado en la Python Console de QGIS se tiene:



No hay comentarios: