[postgis-users] strange behavior in ST_SnaptoGrid

[Mike Taves]

On Mon, 4 Sept 2023 at 00:06, Sajjad Hassany Pazoky
<ijustloveit13 at gmail.com> wrote:
>
> I observed a strange behavior when using ST_SnaptoGrid.
> I don't know whether I am making a mistake or if there is a bug.
> The following code:
> SELECT ST_AsText(ST_SnapToGrid(
>         ST_GeomFromText('LINESTRING(334729.13 4103548.88,334729.12 4103548.53)'),
>         0.001)
>     );
>
> produces:
> LINESTRING(334729.13 4103548.88,334729.12 4103548.5300000003)

This is an annoying quirk of floating point numbers. What is happening
behind the scenes is simple:

value = 4103548.53
grid = 0.001
snapped = round(value / grid) * grid

which results with 4103548.5300000003 which is one floating point unit
too many! There is a nextafter C function that could be used to bring
the result down a notch. Using Python/NumPy the example would be:

np.nextafter(np.rint(value / grid) * grid, -np.inf)

But it's difficult to know when to increment/decrement the float by a
small amount, and which direction. There is a better solution
somewhere, but I haven't figured it out yet.

Note that the actual source code for PostGIS' implementation is:
https://github.com/postgis/postgis/blob/11876ebb5e70dfb3c1b40e48bfc343820e6000d9/liblwgeom/ptarray.c#L2082-L2155