Getting "QgsPolygon: CurvePolygon EMPTY" when extracting the largest part of MultiPolygon using PyQGIS

Question

There is a shapefile with MultiPolygons, see the image below

With this code, I am trying to get the largest part (by area) of each MultiPolygon feature. It utilizes the constGet() method.

layer = iface.activeLayer()

def get_largest_part(feat):
    parts = [part for part in feat.geometry().constGet()]
    largest_part = max(parts, key=lambda x: x.area())
    return largest_part

largest_parts = [get_largest_part(feat) for feat in layer.getFeatures()]

print(largest_parts)

However, the output looks a bit weird. Can somebody explain to me why?

[<QgsPolygon: CurvePolygon EMPTY>, <QgsPolygon: CurvePolygon EMPTY>]

Sometimes when running these both codes my QGIS 3.18.2-Zürich on Windows 8.1 even crashes. The same story happened when using parts() and constGet() on Ubuntu 21.10.

I tried to apply the solution suggested by @KadirŞahbaz in this thread: Getting area of part of MultiPolygon in QGIS, that utilizes the parts() method

layer = iface.activeLayer()

def get_largest_part(feat):
    parts = [part for part in feat.geometry().parts()]
    largest_part = sorted(parts, key=lambda x: x.area())[0]
    return largest_part
    
largest_parts = [get_largest_part(feat) for feat in layer.getFeatures()]
    
print(largest_parts)

however, the output is the same:

[<QgsPolygon: CurvePolygon EMPTY>, <QgsPolygon: CurvePolygon EMPTY>]

Answer 1

I don't know why you get empty polygons, but you can use the following script:

layer = iface.activeLayer()

def get_largest_part(feat):
    geoms = feat.geometry().asGeometryCollection()
    geom = sorted(geoms, key=lambda g: g.area(), reverse=True)[0]
    return QgsGeometry(geom)

largest_parts = [get_largest_part(feat) for feat in layer.getFeatures()]

Answer 2

I've had success by returning a copy of the largest_part geometry within the function.

Unfortunately, I can't really explain why that helps, but I can only assume there is some kind of issue with variable scope.

I guess the QgsAbstractGeometry object gets 'cleared' (to release the resources?) after it leaves the scope of the function, therefore it is an empty polygon if you access it outside of that function. Returning a deep copy with clone() somehow prevents this.

But that's pretty much speculation from my side.

layer = iface.activeLayer()

def get_largest_part(feat):
    parts = [part for part in feat.geometry().constGet()]
    largest_part = max(parts, key=lambda x: x.area())
    return largest_part.clone()

largest_parts = [get_largest_part(feat) for feat in layer.getFeatures()]

print(largest_parts)

# just to highlight geoms on map canvas
geoms = [QgsGeometry(p) for p in largest_parts]
iface.mapCanvas().flashGeometries(geoms, layer.crs())