miércoles, 22 de marzo de 2017

Corrigiendo geometrías con PyQGIS: primera parte

En un post anterior se escribió acerca de la validación y corrección de geometrías y de este link se puede bajar un shapefile de cierta zona de Malasia el cual no permite convertirlo adecuadamente en el polígono esperado a pesar de no presentar, aparentemente, problemas de geometría. Sin embargo, de manera subyacente, presenta muchas inconsistencias.



Así luce al cargarlo a la Map View de QGIS:


y así después de convertirlo a polígono mediante la herramienta "Lines to polygons" de la Processing Tool Box de QGIS:


Si se abre la tabla de atributos de la capa original se observa que es singlepart y los rasgos no se despliegan de manera secuencial; tal como se observa en la imagen siguiente:


El código para validación de geometrías de este post no arroja geometrías inválidas. Por otra parte, en el sitio donde se expuso inicialmente el problema se señala como posible causa la presencia de separaciones entre los respectivos features; uno de alrededor de 20 metros y otro aparentemente milimétrico. Con el código siguiente se van explorar cuantos gaps existen en realidad.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import itertools

layer = iface.activeLayer()

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

n = len(feats)

for i, j in itertools.combinations(range(n), 2):
    distance = feats[i].geometry().distance(feats[j].geometry())
    if distance < 50:
        print i, j, distance

Al ejecutarlo en la Python Console de QGIS se tiene la salida siguiente; donde las dos primeras columnas corresponden a los índices de los rasgos que interactuan y la tercera es la distancia entre ellos:

1
2
3
4
5
6
7
0 3 0.0
0 5 6.32561007326e-08
1 2 0.0
1 4 2.77618203084e-08
2 3 0.0
4 6 3.69805894613e-10
5 6 5.00951291099e-08

La salida sugiere la presencia de 4 gaps muy pequeños pero no hay indicios del que corresponde a alrededor de 20 metros. Si existe, esto señala que uno o varios de los rasgos son a su vez multipartes. Convirtiendo el shapefile original a singlepart con la herramienta correspondiente de la Processing Tool Box, se corrobora que existen rasgos agrupados, 10 en total, como se observa en la imagen a continuación:


Cuando se ejecuta el código anterior para esta capa aparecen nuevos gaps; entre ellos el de 22 metros (para los rasgos 5, 6):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
0 3 0.0
0 7 6.32561007326e-08
1 2 0.0
1 6 2.77618203084e-08
2 3 0.0
4 5 6.09094222774e-08
4 8 3.69805894613e-10
5 6 22.0495367734
7 9 5.00951291099e-08
8 9 9.68909528121e-08

Ahora se van a eliminar los gaps con la herramienta v.clean de la Processing Tool Box con la opcion 'snap' y un umbral de 30 m (que garantiza la elimnación del gap de 22 m); tal como en la imagen siguiente:


Ejecutando nuevamente el código anterior para la capa resultante de v.clean ('4testing_singleparts_clean_snap.shp'), se obtiene:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
0 3 0.0
0 7 0.0
1 2 0.0
1 6 0.0
2 3 0.0
4 5 0.0
4 8 0.0
5 6 0.0
7 9 0.0
8 9 0.0

lo que permite verificar que la herramienta v.clean eliminó totalmente todos los gaps. Sin embargo, esto no era la causa del mal desempeño de "Lines to polygons" porque el resultado es idéntico al aplicarlo a la capa 4testing_singleparts_clean_snap.shp. La resolución de estos problemas los consideraremos en otro post.

No hay comentarios: