5

There is a shapefile with MultiPolygons, see the image below

input

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>]
3
  • 1
    I can't understand what's going on. In Ubuntu 21.10, QGIS always crashes when using parts and constGet. Commented Mar 20, 2022 at 10:15
  • @KadirŞahbaz I will include your comment in my question.
    – Taras
    Commented Mar 20, 2022 at 13:54
  • Crashing occurs usually printing largest_parts in the last line. (rarely in the method's body) QGIS 3.24.0 Commented Mar 20, 2022 at 14:03

2 Answers 2

7

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()]

enter image description here

3
  • Thank you for a good idea. I have a question: have you also tried my code? Did it work?
    – Taras
    Commented Mar 20, 2022 at 7:15
  • It didn't work. I get the same result, too. Commented Mar 20, 2022 at 7:44
  • Can you please tell me, does my code look okay in terms of Python? Is it logical?
    – Taras
    Commented Mar 20, 2022 at 9:32
6

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())

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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