Intersección de rasgos (features) adyacentes de polígonos en PyQGIS

La intersección de rasgos adyacentes en un polígono es posible en PyQGIS a través de los métodos ‘intersects’ e ‘intersection’ de la clase QgsGeometry. El primero permite verificar si existe la intersección y el segundo devuelve un objeto de la clase QgsGeometry que, en el caso de una frontera común entre ambos features, correspondería a una MultiLineString. Existe otro método en QgsGeometry, ‘touches’, que también permite verificar si los rasgos “interaccionan” entre si y se prueba en el código que presento más adelante.

Los objetos del tipo MultiLineString son Multi Part y, por tanto, se representan con un sólo feature en el vectorial; a pesar de que estén constituidos por múltiples segmentos de recta. El vectorial con dos rasgos de la imagen siguiente es el que se utilizó para probar el código.

intersection1

Si tienen más de dos rasgos pueden iterar todas las combinaciones posibles usando la librería itertools para descartar las equivalentes. El código propuesto se encuentra a continuación:

layer = iface.activeLayer()

features = layer.getFeatures()

features = [feature for feature in features]

print "feature 0 touches feature 1?", features[0].geometry().touches(features[1].geometry())

inters = features[0].geometry().intersects(features[1].geometry())

geom1 = features[0].geometry()
geom2 = features[1].geometry()

if inters == True:

    print "Intersection is multipart?", geom1.intersection(geom2).isMultipart()
    
    wkt = geom1.intersection(geom2).exportToWkt()
    
    #creating a memory layer for multi polygon
    crs = layer.crs()
    epsg = crs.postgisSrid()
 
    uri = "MultiLineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
 
    mem_layer = QgsVectorLayer(uri,
                               "multipolyline",
                               "memory")
    
    QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
 
    mem_layer.startEditing()
 
    #Set features
    feature = QgsFeature()
 
    #set geometry
    feature.setGeometry(QgsGeometry.fromWkt(wkt))
    #set attributes values
    feature.setAttributes([1])
    mem_layer.addFeature(feature, True)
 
    #stop editing and save changes
    mem_layer.commitChanges()
    
layer = iface.activeLayer()

feat = layer.getFeatures().next()

print "memory layer is multipart?", feat.geometry().isMultipart()

print wkt

Cuando se ejecuta en la Python Console se produce la memory layer siguiente:

intersection2

donde se verifica por pantalla que es Multi Part. La salida completa se presenta a continuación:

intersection3

Observen que la geometría también se imprime como WKT (Well Known Text) en la Python Console. Si lo desean, pueden prescindir de la memory layer y usar el plugin QuickWKT que tiene el mismo efecto. Solo tienen que hacer copy/paste con la geometría en la ventana del plugin y hacer click en OK. El resultado lo tienen en la imagen siguiente:

intersection4

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

2 respuestas a Intersección de rasgos (features) adyacentes de polígonos en PyQGIS

  1. Pingback: Intersección de polígonos solapados en PyQGIS | El Blog de José Guerrero

  2. Pingback: Plugin de QGIS para determinar fronteras como multipolilíneas | El Blog de José Guerrero

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