Dividir un polígono por otro mediante PyQGIS

Supongamos que se tiene una capa polígono que representa las unidades administrativas de una región y sobre ella colocamos otra que corresponde a lotes de terreno. Es probable que algunas de estas últimas abarquen más de una zona de la primera, tal como se presenta en la imagen siguiente, y por tanto nos interese reflejar esa división.

dividir1

La herramienta de geoproceso que nos permite dar cuenta de esa situación es la intersección y el código a continuación nos sirve para ese propósito:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feats1 = [ feat for feat in layers[0].getFeatures() ]

feats2 = [ feat for feat in layers[1].getFeatures() ]

new_features = []

for feats in feats1:
    geom1 = feats.geometry()
    for feats in feats2:
        geom2 = feats.geometry()
        if geom1.intersects(geom2):
            geom3 = geom1.intersection(geom2).exportToWkt()
            new_feat = QgsFeature()
            new_feat.setGeometry(QgsGeometry.fromWkt(geom3))
            new_features.append(new_feat)

#creating a memory layer for polygons
crs = layers[0].crs()
epsg = crs.postgisSrid()
      
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
      
mem_layer = QgsVectorLayer(uri,
                           "split_polygon",
                           "memory")
      
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

mem_layer.startEditing()

n = len(new_features)

for i in range(n): 
    #set attributes values
    new_features[i].setAttributes([i])
    mem_layer.addFeature(new_features[i], True)

#stop editing and save changes
mem_layer.commitChanges()

Después de ejcutado el código anterior en la Python Console de QGIS se obtiene lo siguiente:

dividir2

Se han seleccionado arbitrariamente algunos de los rasgos objeto de división que corroboran que el código funcionó de la manera esperada.

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

3 respuestas a Dividir un polígono por otro mediante PyQGIS

  1. Pingback: Seleccionar rasgos adyacentes en un polígono mediante PyQGIS | El Blog de José Guerrero

  2. Giannina B dijo:

    Gracias por el artículo José. Qué pasa si añado un polígono desde Google Earth con extensión .kml y lo paso a .shp en qgis… es posible ensamblar a este nuevo polígono con las capas ya existentes?

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