2

I'm trying to intersect two polygons with intersection() function of QgsGeometry (pyqgis-3). The problem is that the result is a multipart polygon instead a single part polygon and QGIS doesn't allow me to save the ouput layer. The output layer needs to be a Polygon, not Multipolygon. I need to work at geometry level, I followed this topic but the problem already exist.

This is my python code:

l1_list = QgsProject.instance().mapLayersByName('Buffer_10km')
l2_list = QgsProject.instance().mapLayersByName('municipality')

l1 = l1_list[0]
l2 = l2_list[0]

prv1 = l1.dataProvider()
prv2 = l2.dataProvider()

# Create a memory layer to store the result
outlayer = QgsVectorLayer('Polygon?crs=' + prv1.crs().toWkt(), "result", "memory")
outpr = outlayer.dataProvider()

feats_l1 = l1.getFeatures()
list_geom_l1 = []
for feat_l1 in feats_l1:
    print('id feature: ' + str(feat_l1.id()))
    g_l1 = feat_l1.geometry()
    list_geom_l1.append(g_l1)

outlayer.startEditing()
outlayer. dataProvider().addAttributes([QgsField('code', QVariant.String)])

for feature in prv2.getFeatures():
    g = feature.geometry()
    g_tmp = g.intersection(list_geom_l1[0])
    print('Multipart?: ' + str(g_tmp.isMultipart()))
    new_feat = QgsFeature()
    new_feat.setGeometry(g_tmp)
    outlayer.addFeature(new_feat)

outlayer.commitChanges()
outlayer.updateExtents()
QgsProject.instance().addMapLayer(outlayer)

The reference system of input layers are the same (WGS84). Below there is an image of the layers.

enter image description here

  • 2
    The intersection is a multipart polygon with two parts. Do you want turn it into two single part polygons or only keep one part? – BERA Aug 5 at 10:47
  • Turn it into two single part polygons – PaolaB Aug 5 at 10:49
1

I found a solution thanks to osgeo library that allows to interact with the single parts of a multipolygon geometry. Help me the answer of Barret in this topic.

from osgeo import ogr
...
for feature in prv2.getFeatures():
    g = feature.geometry()
    #Osgeo library
    g1wkt = ogr.CreateGeometryFromWkt(g.asWkt())
    g2wkt = ogr.CreateGeometryFromWkt(list_geom_l1[0].asWkt())
    g_tmp = g1wkt.Intersection(g2wkt)    
    if g_tmp.GetGeometryName() == 'MULTIPOLYGON':
        for geom_part in g_tmp:
            gw = geom_part.ExportToWkt()
            gg = QgsGeometry.fromWkt(gw)
            #go back to pyqgs geometry
            new_feat = QgsFeature()
            new_feat.setGeometry(gg)
            outlayer.addFeature(new_feat)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.