Multipart geometry after PyQGIS-3 intersection()

Question

I'm trying to intersect two polygons with intersection() function of QgsGeometry. The problem is that the result is a multipart polygon instead a single part polygon and QGIS doesn't allow me to save the output 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.

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 05 '19 at 10:47

PaolaB · Answer 1 · 2019-08-12T14:41:32.597

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)

Multipart geometry after PyQGIS-3 intersection()

1 Answers1