Reducing precision results in empty polygon unnecessarily (GEOSGeom_setPrecision)

[jorisvandenbossche]

Original report: https://github.com/shapely/shapely/issues/2025

Reproducer with geosop using latest main:

$ geosop -a "POLYGON ((1 0, 0 6, 1 3, 1 0))" reducePrecision 1.0
POLYGON EMPTY

In this case, it seems the precision of 1.0 should be able to perfectly preserve the input geometry, however an empty polygon is returned. Is there an explanation for why this would be the expected result, or is that a bug?

The KeepCollapsed version also results in an empty polygon, but Pointwise preserves the input:

$ geosop -a "POLYGON ((1 0, 0 6, 1 3, 1 0))" reducePrecisionKeepCollapsed 1.0
POLYGON EMPTY
$ geosop -a "POLYGON ((1 0, 0 6, 1 3, 1 0))" reducePrecisionPointwise 1.0
POLYGON ((1 0, 0 6, 1 3, 1 0))

[dr-jts]

This is because the algorithm used to reduce precision actually snap-rounds all vertices and edges to a grid with the specified precision. In this case that causes the longer line to snap to the nearby vertex, and thus effectively collapse the polygon.

I agree this behaviour seems a bit counter-intuitive. I'm not sure it's possible to change it, at least not using the current approach. And in fact it might be that trying to come up with an approach to change this result might produce unwanted results in other situations.

[jorisvandenbossche]

In this case that causes the longer line to snap to the nearby vertex, and thus effectively collapse the polygon.

So it's the line that snaps? Because the vertices itself, 0 6, fall on the grid and wouldn't snap? I would naively (not having actual knowledge about the details of the algorithm) assume that only the vertices are snap rounded, but that's a wrong mental model?

[dr-jts]

So it's the line that snaps? Because the vertices itself, 0 6, fall on the grid and wouldn't snap?

Correct. Vertices are rounded to the nearest grid point, and lines are snapped to vertices which lie in grid cells that the line intersects. This is done to avoid situations where a vertex moves and crosses a line.

For example, in the following case the vertices lie on grid points, but the longest line is snapped to the nearby vertex and causes the polygon to be split in two:

bin/geosop -a "POLYGON ((9 1, 1 5, 1 1, 4 3, 4 1, 9 1))" reducePrecisionKeepCollapsed 1.0

MULTIPOLYGON (((9 1, 4 1, 4 3, 9 1)), ((1 1, 1 5, 4 3, 1 1)))

I suppose it might be possible to snap lines only if they are crossed by vertices that have moved. Have to think about that.

[jorisvandenbossche]

Got it. Thanks for the explanation!

[theroggy]

If you would change this behaviour, please make it configurable, as I use this a lot to clean up slivers and almost dangling nodes after applying overlay operations...

[dr-jts]

If you would change this behaviour, please make it configurable, as I use this a lot to clean up slivers and almost dangling nodes after applying overlay operations...

Good to know there is a clear use case for this behaviour. No plans to change anything at this time, but if we do then we will provide an option.

[bretttully]

Thanks for the explanation @dr-jts -- this result is counter intuitive to me (hence the original shapely issue), but we also share the use case mentioned by @theroggy so it might actually be a good thing now we know how it works!

Is there a reference you can point to for the algorithm that is implemented in GEOS/JTS? Or is this one of your own creation? I'd like to go deeper and try and build out some visual documentation.

[bretttully]

Looking at the cgal docs and they distinguish between snap rounding and iterated snap rounding with a nice picture to show the difference. https://doc.cgal.org/latest/Snap_rounding_2/index.html

Am I correct in understanding your comment above that GEOS is implementing ISR so any edge that crosses a hot pixel has a node added at that pixel? That would explain why as the polygon gets a bigger aspect ratio it collapses to a single line/empty polygon.

[dr-jts]

Looking at the cgal docs and they distinguish between snap rounding and iterated snap rounding with a nice picture to show the difference. https://doc.cgal.org/latest/Snap_rounding_2/index.html

That's a good reference. The original papers are ok too (although they leave out a few small but key details).

Am I correct in understanding your comment above that GEOS is implementing ISR so any edge that crosses a hot pixel has a node added at that pixel? That would explain why as the polygon gets a bigger aspect ratio it collapses to a single line/empty polygon.

No, there is no iteration in the GEOS algorithm. It uses standard snap-rounding.

[dr-jts]

Is there a reference you can point to for the algorithm that is implemented in GEOS/JTS? Or is this one of your own creation? I'd like to go deeper and try and build out some visual documentation.

The concept of snap-rounding is explained fairly well in the CGAL docs mentioned above: https://doc.cgal.org/latest/Snap_rounding_2/index.html

The original papers are available online as well, e.g.

• M. Goodrich, L. J. Guibas, J. Hershberger, and P. Tanenbaum. Snap rounding line segments efficiently in two and three dimensions

There's still quite a bit of details that need to be added to get a fully functioning, performant algorithm. The only reference for those is the code.

[bretttully]

This makes a lot of sense now, thank you!

For those who see this in the future, the following plot helped me in this particular case. You can see in the latter two columns, the hypotenuse intersects with a hot pixel -- snap rounding then turns the hypotenuse into a polyline with a vertex at that hot pixel, and thus the polygon collapses to zero area.

import matplotlib.pyplot as plt
import numpy as np
import shapely

fig, axes = plt.subplots(ncols=4, figsize=(10, 20))
for i, ax in zip(range(3, 7), axes, strict=True):
    background = np.zeros((7, 3), dtype=np.uint8)
    p = shapely.Polygon([(1,0), (0, i), (1, 3), (1,0)])
    for x, y in p.exterior.coords:
        background[int(y), int(x)] = 255
    ax.imshow(background, cmap='gray')
    ax.plot(*p.exterior.xy, color='red', linewidth=2, marker='o')
    ax.set_aspect('equal')
    ax.set_title(f'Collapses? {shapely.set_precision(p, 1).is_empty}')

plt.show()
plt.close(fig)

[theroggy]

Documentation wise, for a while now there has been a PR in the works to document (amongst others) this behaviour in the shapely docs: https://github.com/shapely/shapely/pull/1964/files

Possibly the info in this thread can offer info to further finetune this PR. Not sure which level of detail is expected/wanted in the shapely docs.