Intersecciones con shapely en shapefiles con número de features elevado

Una de las desventajas que tienen los lenguajes interpretados es la velocidad de proceso y eso de pone de manifiesto cuando nos toca hacerlo con shapefiles con gran número de features; tal como en el del ejemplo siguiente. El shapefile tiene 20.730 features y el número de combinaciones, tomadas de 2 en 2 para evaluar las primeras intersecciones, es de 214.856.085.

Por esta razon, no es extraño la lentitud de la herramienta optimizada de la ‘processing toolbox’ para obtener las ‘primeras intersecciones’; sin pensar en la aún mayor lentitud de un script propio usando los métodos de QgsGeometry (‘intersects’ e ‘intersection’).

Sin embargo, el uso de los métodos ‘unary_union’ y ‘polygonize’ del módulo Python de shapely resultan muy promisorios en la obtención de todas las intersecciones (incluyendo las intersecciones de intersecciones); tal como se presenta en el script siguiente:

import fiona
from shapely.geometry import shape, LineString
from shapely.ops import unary_union, polygonize
import time

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

collection = []

start = time.time()

print "wait..."

for i,item in enumerate(c):
    geom = shape(item['geometry'])
    collection.append(geom)

rings = [ LineString(pol.exterior.coords) for pol in collection ]

union = unary_union(rings)

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

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

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

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

prov = mem_layer.dataProvider()

for i,feat in enumerate(new_intersections):
    feat = QgsFeature()
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(new_intersections[i]))
    prov.addFeatures([feat])

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

end = time.time()

time_tot = end - start

print "total time = %.4f s" % time_tot

Después de alrededor de unos 18 minutos de proceso se obtiene una memory layer con 57.401 features; correspondientes a todas las intersecciones.

En la imagen siguiente se ha seleccionado arbitrariamente un rasgo en la capa original con 50 % de transparencia. Esto permite detectar la superposición de otros features susceptibles de intersectar con el ya seleccionado.

intersection1

En la imagen siguiente se visualiza la memory layer resultante, también con 50 % de transparencia, donde se permite apreciar que el rasgo seleccionado en el shapefile original se ha convertido en un “mosaico” de 7 features correspondiente a intersecciones de intersecciones.

intersection2

Anuncios
Esta entrada fue publicada en PyQGIS, 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