Rings order in MultiPolygon WKT string may leave overlapping parts when deleting holes (native:deleteholes)

jfbourdon commented 3 years ago

What is the bug or the crash?

I'm not so sure that this is a bug or a feature request but the issue it that when manually digitalizing a MultiPolygon, multiple level of nested rings don't seem to be store correctly or at least not in a way that some processing tools are expecting.

The following MultiPolygon, when manually drawn is defined as: MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2))) multipart

When trying to fill the hole using native:deleteholes, the hole is filled, but the island is not dissolved, resulting in an invalid feature (overlapping parts of the same feature): MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2))) multipart_notdissolved

Instead of the expected: MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0)),((6 2, 6 3, 7 3, 7 2, 6 2))) multipart_dissolved

However, if the WKT string is defined as 'MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1),(2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))' then the island is dissolved.

Steps to reproduce the issue

MultiPolygon with overlapping island after native:deleteholes

vlayer = QgsVectorLayer("multipolygon?crs=EPSG:3857", "test", "memory")
geom = QgsGeometry().fromWkt('MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))')
fet = QgsFeature()
fet.setGeometry(geom)
vlayer.dataProvider().addFeature(fet)

params = {'INPUT':vlayer, 'MIN_AREA':10, 'OUTPUT':'TEMPORARY_OUTPUT'}
vlayer_filled = processing.run("native:deleteholes", params)["OUTPUT"]

vlayer_filled.getGeometry(1).asWkt()
# 'MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))'

MultiPolygon with dissolved island after native:deleteholes

vlayer = QgsVectorLayer("multipolygon?crs=EPSG:3857", "test", "memory")
geom = QgsGeometry().fromWkt('MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1),(2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))')
fet = QgsFeature()
fet.setGeometry(geom)
vlayer.dataProvider().addFeature(fet)

params = {'INPUT':vlayer, 'MIN_AREA':10, 'OUTPUT':'TEMPORARY_OUTPUT'}
vlayer_filled = processing.run("native:deleteholes", params)["OUTPUT"]

vlayer_filled.getGeometry(1).asWkt()
# 'MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0)),((6 2, 6 3, 7 3, 7 2, 6 2)))'

Versions

QGIS version 3.20.0-Odense QGIS code revision decaadbb31 Qt version 5.15.2 Python version 3.9.5 GDAL/OGR version 3.3.0 GEOS version 3.9.1-CAPI-1.14.2 Windows 10 Version 1809

Additional context

I'm not sure of the fix:

That native:deleteholes should automatically remove overlapping parts
That nested ring should be written inside the same parenthesis even though the OCG spec don't require that (as far as I understand it)
```
MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1)),((2 2, 2 3, 3 3, 3 2, 2 2)))
MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1),(2 2, 2 3, 3 3, 3 2, 2 2)))
```
Interesting read on GIS StackExchange on the matter: Order of polygon vertices in general GIS.

jfbourdon commented 3 years ago

I should also add that the order of rings in the WKT has an impact in the ability of native:deleteholes to account for the area occupied the island inside the ring. In fact, it seems that the area of islands is never considered.

The following MultiPolygon has a inner ring of 9 unit² with an island of 1 unit² resulting in a hole of 8 unit². So if I set the MIN_AREA parameter to 8.5, I expect the hole to disappear but it won't.

vlayer = QgsVectorLayer("multipolygon?crs=EPSG:3857", "test", "memory")
geom = QgsGeometry().fromWkt('MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))')
fet = QgsFeature()
fet.setGeometry(geom)
vlayer.dataProvider().addFeature(fet)

params = {'INPUT':vlayer, 'MIN_AREA':8.5, 'OUTPUT':'TEMPORARY_OUTPUT'}
vlayer_filled = processing.run("native:deleteholes", params)["OUTPUT"]

vlayer_filled.getGeometry(1).asWkt()
# 'MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))'

Interestingly, if I put all the nested ring inside the same parenthesis as I suggest, the hole also stays but the island is removed.

vlayer = QgsVectorLayer("multipolygon?crs=EPSG:3857", "test", "memory")
geom = QgsGeometry().fromWkt('MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1),(2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))')
fet = QgsFeature()
fet.setGeometry(geom)
vlayer.dataProvider().addFeature(fet)

params = {'INPUT':vlayer, 'MIN_AREA':8.5, 'OUTPUT':'TEMPORARY_OUTPUT'}
vlayer_filled = processing.run("native:deleteholes", params)["OUTPUT"]

vlayer_filled.getGeometry(1).asWkt()
# 'MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1)),((6 2, 6 3, 7 3, 7 2, 6 2)))'

As a side note, I find misleading the parameter MIN_AREA as it rather define the maximum area of holes to be deleted (the value itself being excluded).

lbartoletti commented 3 years ago

Thank you for this very clear and detailed report, it is much appreciated.

I am not sure if this can be considered a bug. Indeed, the algorithm does remove the hole, that's what it is asked. However, you're right it produces an invalid geometry. There are some others algorithm which can produce invalid geometry (I think about snap to grid, subdivide, etc)

This case is simple to solve, since we can extract and merge polygon. Or use the magic buffer(0) trick, example:

geom = QgsGeometry().fromWkt('MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0),(1 1, 4 1, 4 4, 1 4, 1 1)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))')

# removeInteriorRings result:
geom.removeInteriorRings()
# <QgsGeometry: MultiPolygon (((0 0, 0 5, 5 5, 5 0, 0 0)),((2 2, 2 3, 3 3, 3 2, 2 2)),((6 2, 6 3, 7 3, 7 2, 6 2)))>

# valideGeometry result
geom.removeInteriorRings().validateGeometry()
[<QgsGeometry.Error:  Le polygone 1 est à l'intérieur du polygone 0>]

# using unaryUnion:
QgsGeometry().unaryUnion([geom.removeInteriorRings()])
<QgsGeometry: MultiPolygon (((0 5, 5 5, 5 0, 0 0, 0 5)),((7 3, 7 2, 6 2, 6 3, 7 3)))>

# using buffer(0)
geom.removeInteriorRings().buffer(0, 0)
# <QgsGeometry: MultiPolygon (((6 2, 6 3, 7 3, 7 2, 6 2)),((0 0, 0 5, 5 5, 5 0, 0 0)))>

# I noticed that makeValid doesn't work for this case
geom.removeInteriorRings().makeValid()
#<QgsGeometry: MultiPolygon (((5 5, 5 0, 0 0, 0 5, 5 5),(3 3, 2 3, 2 2, 3 2, 3 3)),((7 3, 7 2, 6 2, 6 3, 7 3)))>

qgis / QGIS