Open theroggy opened 2 weeks ago
This is an issue in JTS as well. It's a robustness issue in the overlay code.
No fix at the moment. It does work using Snap Rounding, with a fairly large scale factor (10^7). Not sure if that's available in Shapely.
It does work using Snap Rounding, with a fairly large scale factor (10^7). Not sure if that's available in Shapely.
I don't immediately find anything regarding snaprounding in the shapely API, nor in the GEOS C API, so I suppose not.
I have 4 more cases that also lead to "TopologyException: Ring edge missing at...", and 3 that give a "TopologyException: found non-noded intersection between ..." error.
Is it useful to find the exact culprits for those cases as well (no problem to do so, it takes about 5 to 10 minutes per case to seperate them from the 3 million input polygons in total) or is it very likely that they all lead back to the same issue?
It does work using Snap Rounding, with a fairly large scale factor (10^7). Not sure if that's available in Shapely.
I don't immediately find anything regarding snaprounding in the shapely API, nor in the GEOS C API, so I suppose not.
It looks like the Shapely overlay functions accept a gridSize
parameter.
Is it useful to find the exact culprits for those cases as well (no problem to do so, it takes about 5 to 10 minutes per case to seperate them from the 3 million input polygons in total) or is it very likely that they all lead back to the same issue?
It would be interesting to see some other failure cases. They may have different causes, which may (or may not) be easier to fix.
OK, I further automated the discovery of the extra cases... each .csv contains some WKTs that when unioned trigger an error.
error_6.csv error_0.csv error_1.csv error_2.csv error_3.csv error_4.csv error_5.csv
It looks like the Shapely overlay functions accept a
gridSize
parameter.
Ah, oops, ok... grid_size is synonymous for snap-rounding... I should have known, sounds logical (blushing).
OK, I further automated the discovery of the extra cases... each .csv contains some WKTs that when unioned trigger an error.
Thanks. I confirm I see the errors as well. I can't offer any fast resolution for this, unfortunately. But you might try either using snap-rounding (via specifying a grid size in Shapely), or doing some amount of precision reduction on the coordinate values (which will likely move coordinates enough to avoid these problems.
Since you're using Shapely, can you catch the error and use snap-rounding only on the cases which cause the error? (It may be that this strategy can be built into the Unary Union code itself - but I haven't figured out yet how to do this in a way that won't reduce performance across the board.)
I was actually allready applying gridsize=0.01
(the data is projected, so 0.01 meter is still plenty precise), but I'm kind of using geopandas, and in geopandas the gridsize
parameter isn't exposed yet... so I was just using set_precision
after doing the unary_union
... which obviously doesn't help to avoid this error.
Because of your advice that using snap-rounding/gridsize in the unary_union
can avoid the error, I now made a version where the gridsize
parameter does trickle down so it is directly applied in the unary_union
, and this solves all of the above issues.
Thanks a lot already!
Because of your advice that using snap-rounding/gridsize in the
unary_union
can avoid the error, I now made a version where thegridsize
parameter does trickle down so it is directly applied in theunary_union
, and this solves all of the above issues.Thanks a lot already!
Nice! Glad that helped.
I'll keep thinking about a more fundamental fix that will solve this problem for all users of JTS and GEOS.
The fundamental problem is that it can happen that unioning two geometries produces a result which has invalid polygonal elements. This probably happens due to almost coincident linework causing the overlay topology determination to be incorrect.
Here is a small reproducer example:
POLYGON ((39893.875 197462.125, 39900 197400, 39820.57741178792 197500, 39848.89740965763 197500, 39893.875 197462.125), (39891.625 197448.375, 39872.007923812 197457.49692859003, 39862.259636527146 197475.4252251145, 39895.67580000311 197413.96880000085, 39891.625 197448.375))
POLYGON ((39866.13 197473.47, 39868.73410547314 197475.00383519547, 39866.12902312726 197473.4659950994, 39863.954506024966 197472.3081521314, 39863.95342412875 197472.3101418713, 39866.13 197473.47))
Input, showing very narrow B polygon touching a gore in A
Magnified Topology of input
Currently the union of A and B produces:
GEOMETRYCOLLECTION (POLYGON ((39900 197400, 39820.57741178792 197500, 39848.89740965763 197500, 39893.875 197462.125, 39900 197400), (39872.007923812 197457.49692859003, 39895.67580000311 197413.96880000085, 39891.625 197448.375, 39872.007923812 197457.49692859003)), POLYGON ((39863.95356604543 197472.3098808691, 39863.95342412875 197472.3101418713, 39866.13 197473.47, 39868.73410547314 197475.00383519547, 39866.12902312726 197473.4659950994, 39863.954506024966 197472.3081521314, 39863.95356604543 197472.3098808691)), LINESTRING (39863.95342412875 197472.3101418713, 39862.259636527146 197475.4252251145))
Result containing a nested polygon element
The polygonal elements in this geometry overlap, so their combination is an invalid MultiPolygon. CascadedUnion
then unions that with further geometries. This produces a TopologyException
due to the invalid input.
A fix is to have CascadedUnion
ensure that polygonal inputs extracted from mixed GeometryCollection
s are valid, and if not to run GeometryFixer.fix
on them. This might not cause too much performance degradation, since mixed-dimension GCs should be relatively rare.
This is similar to #942. And like that situation, this can't be solved by simple area or envelope check heuristics (since the results lie within the inputs).
Doing a union of a combination of 4 specific geometries leads to a TopologyException: Ring edge missing at ... .
They aren't the prettiest geometries, and they overlap... but they are all valid, so...
Script to reproduce (using shapely 2.0.6 with geos 3.13)