Opposite intersection and difference results

[mortewle]

This might be a difficult-to-solve issue, but I get some weird results from polygon overlays. One example follows.

I have two valid polygons (attatched), where geom2 has a very thin hole in it. The intersection with geom1 should be everything except the hole, but the resulting geometry is only the hole. The opposite is the case after difference overlay.

I haven't been able to install geosop, but have tested that this happens in QGIS with GEOS 3.12, as well as in QGIS with GEOS 3.11.2, shapely with GEOS 3.11.1 and sf with GEOS 3.7.

geom1.txt

geom2.txt

[pramsey]

So, that hole is about a nanometer wide. This kind of argues for ST_MakeClean, again, something to take "valid" geometries and remove ridiculous components like "holes that are less than a molecule wide", and various other structurally useless things like spikes and gores and ultra-fine zigzags. We have some prior art of these kinds of inputs, but dealing with them is still an unsolved problem. You might find that using the fixed precision variant of ST_Intersection results in a "more correct" answer (hopefully) although at the cost of very slightly rounding all your vertices in the process. If you don't have exact edge-matching in your data set it might be a useful approach.

[mortewle]

Thanks for the answer @pramsey!

Context: we are in the process of refactoring a statistics production previously done without trouble in ArcGIS, so these kinds of differences are interesting.

I realize this is a tiny hole and easily solvable. This is one of many examples.

The more problematic issues (should maybe be a separate issue?) are when it’s spikes, as you mention. The best solution I have, small negative buffer, often enough raises a topology exception. I would happily take some tips on how to deal with this.

[pramsey]

In this case as well, I'd be interested to know if something like GEOSIntersectionPrec() helps. If you're running in PostGIS, the three-parameter version of ST_Intersection(). If the tolerance is big enough it should collapse some of these things down.

[mortewle]

Not running PostGIS, but geopandas/shapely. The shapely overlay operations offer a grid_size parameter, which looks like what GEOSIntersectionPrec is doing.

I have tried this previously, with no success (a lot of geometries disappeared with no warning), but I will try again.

[pramsey]

Hum, disappearing geometries seems bad, particularly since the *Prec code path is supposed to drop us into our most reliable algorithms (with the caveat that every vertex ends up snap rounded). I wonder if it's because after rounding the geometry is invalid (self-intersecting ring). @dr-jts ?

[dr-jts]

I'll have to dig into this to see what's going on. This feels like a known issue where numerical robustness failures cause the topology to be computed incorrectly, and the safeguard heuristics in OverlayNG don't detect the problem.

[dr-jts]

@mortewle Do you have some examples of "spike" cases that are causing problems?

[dr-jts]

In this case as well, I'd be interested to know if something like GEOSIntersectionPrec() helps.

I tried using JTS OverlayNG with fixed-precision Snap-Rounding (scale factor = 10^6). This causes the overlay to work correctly. (I'll confirm this in GEOS as well).

So this is a potential fix, although not ideal, since Snap-Rounding causes all coordinates to be rounded. And to use it to resolve occasional errors requires that it be used for all cases.

[mortewle]

Yes, this example works correctly with 1e-6 rounding for me as well, thanks! The rounding of the coordinates is not a big issue in this production, but it might be in other cases.

The program takes some time to complete, but I am yet to encounter any dissappearing or doublet polygons that are wider than the grid size. There might have been some other issue when I last tested this approach.