Usando Fiona y Shapely en la Python Console de QGIS (GNU/Linux Debian)

Hace unos días publiqué un primer post acerca de ciertos módulos complementarios (GeoPandas) a la labor que puede realizarse con PyQGIS. Hoy me tocará escribir sobre Fiona y Shapely.

Fiona es un contenedor Python para funciones de acceso a datos vectoriales de la biblioteca OGR. Los registros de datos de archivos vectoriales los lee en formato GeoJSON y contiene métodos que permite hacerlos disponibles para Shapely (o PyQGIS) o escribir el mismo tipo de asignaciones como registros en archivos nuevos. Por otra parte, Shapely es un paquete de Python para el análisis de teoría de conjuntos y manipulación de las características planares de geometrías mediante funciones de la bien conocida librería GEOS.

Para probar las librerías voy a usar el shapefile de la imagen siguiente pero, en realidad, no necesito tener cargada la capa en el Map Canvas.

fiona1

Las siguientes líneas de código ejemplifican la lectura del shape desde su ruta en disco y el uso del método ‘shape’ de ‘shapely.geometry’ para hacer disponibles los rasgos como geometrías listas para su manipulación. El loop permite crear las listas ‘points’ y ‘collection’.

import fiona
from shapely.geometry import shape

c = fiona.open('/home/zeito/pyqgis_data/circles.shp')

points = []
collection = []

for item in c:
    geom = shape(item['geometry'])
    points.append(geom.centroid.wkt)
    collection.append(geom)

Si se imprimen ambas listas en la Python Console de QGIS se obtiene:

>>>points
['POINT (383088.9326728958 4438831.545159166)', 
'POINT (387356.386428252 4465503.131130142)', 
'POINT (410614.0093949448 4444165.862353361)']
>>>collection
[<shapely.geometry.polygon.Polygon object at 0x8bf7c30c>, 
<shapely.geometry.polygon.Polygon object at 0x8bf7ccec>, 
<shapely.geometry.polygon.Polygon object at 0x8bf7c84c>]

lo que permite corroborar que cada elemento de points está en el formato WKT (Well Known Text) y los de collection son objetos de shapely.geometry.polygon.Polygon.

Si quiero crear una memory layer con PyQGIS requiero la proyección y esta debe ser obtenida con Fiona; Shapely sólo maneja geometrías (nada de proyecciones y sus transformaciones). Por tanto, la memory layer de centroides almacenados en la lista ‘points’ se obtiene con:

epsg = int(c.crs.values()[0].split(':')[1])

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

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

prov = mem_layer.dataProvider()

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

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

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

que se visualiza de esta manera en el Map Canvas:

shapely1

Importando el módulo itertools es posible obtener todas las posibles intersecciones en el mismo shapefile sin repetición; tal como se encuentra a continuación.

import itertools
intersections = [poly1.intersection(poly2).wkt 
                 for poly1,poly2 in  itertools.combinations(collection, 2) 
                 if poly1.intersects(poly2)]

La lista ‘intersections’ contiene todos los WKT correspondientes a esas geometrías que resultaron de la intersección de aquellas almacenadas en la lista ‘collection’. El despliegue en la próxima imagen es posible gracias al plugin QuickWKT.

shapely2

No obstante, se observa que las intersecciones resultantes no ponen de manifiesto las posibles intersecciones entre ellas mismas. Para lograr ésto Shapely pone a la disposición el método ‘unary_union’ de shapely.ops el cual permite “cortar” las LinearRings/LineString de cada intersección las cuales, en el caso de las curvas cerradas, pueden ser poligonizadas con el método ‘poligonize’ del mismo módulo. En este caso:

from shapely.ops import unary_union, polygonize
rings = [ LineString(pol.exterior.coords) for pol in collection ]

union = unary_union(rings)

new_intersections = [geom.wkt for geom in polygonize(union)]

lo cual se visualiza rápidamente con el plugin QuickWKT de QGIS de la siguiente manera:

shapely3

Esta entrada fue publicada en Código Python, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s