Calculate difference with GeometryCollection containing overlapping polygons

[theroggy]

I noticed that calculating the difference with a GeometryCollection that contains overlapping polygons results in the following error: TopologyException: side location conflict at 5 0. This can occur if the input geometry is invalid.

I found some threads from a few years back that stated already that this was not supported at that moment and that it might be considered in the future. I wonder, are there plans to support this in the future?

Script to reproduce in python:

import shapely

# Difference of geom: 1 polygon, geom: GeometryCollection of 2 overlapping polygons
# Result: TopologyException
geom = shapely.Polygon([(0, 0), (50, 0), (50, 50), (0, 50), (0, 0)])
geom_collection = shapely.GeometryCollection(
    [
        shapely.Polygon([(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]),
        shapely.Polygon([(5, 0), (20, 0), (20, 10), (5, 10), (5, 0)]),
    ]
)
print(f"versions: shapely: {shapely.__version__}, geos-C-api: {shapely.geos_capi_version_string}")

print(f"is simple poly valid? {shapely.is_valid_reason(geom)}")
print(f"is collection valid? {shapely.is_valid_reason(geom_collection)}")

try:
    result = shapely.difference(geom, geom_collection)
except Exception as ex:
    print(f"Exception raised: {ex}")

Output:

versions: shapely: 2.0.1, geos-C-api: 3.12.0-CAPI-1.18.0
is simple poly valid? Valid Geometry
is collection valid? Valid Geometry
Exception raised: TopologyException: side location conflict at 5 10. This can occur if the input geometry is invalid.

[pramsey]

What GEOS version? Real support of GeometryCollection overlay only landed in 3.12

[theroggy]

What GEOS version?

Apparently I first tested with 3.11.2 but now I also tested with 3.12.0 which gives the same result. I added the logging I used to verify the version in the script above for reference.

Real support of GeometryCollection overlay only landed in 3.12

Ah, great that the effort to support it has landed already!

[theroggy]

I did some further digging, and found the following relevant changes in 3.12 in the changelog:

• #716: here the documentation was also updated making explicid what input is supported now in OverlayNG: OverlayNG class reference Based on that, overlapping polygons in GeometryCollections aren't supported (yet).
• #797: if I understand correctly this is a change to RobustNG but it's not entirely clear to me what impact it had on the above case.

Because I'm using geos via shapely, it is actually the C API that is being called, so GEOSDifference_r or GEOSDifferencePrec. I'm not sure if this leads to only OverlayNG being called and/or if it has an (automatic) fallback to RubustNG.

[pramsey]

You'll get slightly different code paths depending on if you use GEOSDifference or GEOSDifferencePrec. The former will go via Geometry::difference which picks up the handling of mixed collections built into HeuristicOverlay. The latter will go directly into OverlayNGRobust::Overlay, which does not have mixed collection handling.

This is clearly SubOptimal(tm) and perhaps in 3.13 we should move the mixed dimension handling into OverlayNG so that all ops pick up that logic. cc/ @dr-jts