On Mon, 4 Sept 2023 at 00:06, Sajjad Hassany Pazoky
<ijustlovei...@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
_______________________________________________
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users
