Intersección de líneas con PyQGIS

Si observan las opciones de geoproceso de QGIS (Vector -> Geoprocessing Tools) intuirán, por los iconos, que éstas están referidas a vectoriales tipo polígono. No obstante, es obvio vislumbrar que la intersección de, por ejemplo, dos rectas no paralelas tendrán como intersección un punto en el plano. Por tanto, en los SIG podemos esperar que polilíneas cuyos segmentos de recta cambian de dirección varias veces o segmentos de recta suficientemente largos puedan tener más de un punto de intersección con otras polilíneas o, incluso, geometrías de tipo polígono.

Para solventar esta falta de aplicaciones relacionadas con geoproceso podemos apelar nuevamente a los metodos ‘intersects’ e ‘intersection’ de QgsGeometry en PyQGIS. Para probar mi procedimiento creé el vectorial de líneas de la imagen siguiente que contiene 10 rasgos (features). Es apropiado para utilizar la librería itertools y así encontrar de manera eficiente (sin repetición) todas las combinaciones 2 a 2 que producen puntos en su intersección.

intersection_lines1

El código propuesto es el siguiente:

from utils import createMemoryLayer #From my own library

import itertools

layer = iface.activeLayer()

features = [ feature for feature in layer.getFeatures() ]

n = layer.featureCount()

list = range(n)

points = [features[i].geometry().intersection(features[j].geometry()).asPoint()
          for i,j in itertools.combinations(list, 2)
          if features[i].geometry().intersects(features[j].geometry())]

crs = layer.crs()
epsg = crs.postgisSrid()
 
createMemoryLayer(QGis.Point, epsg, points)

Utiliza la función ‘createMemoryLayer’, que está dentro de mi propio módulo de funciones que he dado en llamar utils, y que es válida de usar siempre que las listas de QgsPoint que se pasan como parámetro se ajusten a las geometrías de tipo punto, línea o poligono (no como WKT).

Después de ejecutado el código anterior en la Python Console, se obtiene la memory layer de la imagen siguiente que contiene los 24 puntos esperados para la intersección de los 10 features.

intersection_lines2

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