Unexpected rounding from ST_SnapToGrid

[michaudm]

I would expect that ST_Equals(ST_GeomFromText('POINT(311.4 0)'), ST_SnapToGrid(ST_GeomFromText('POINT(311.4 0)'), 0.1));

but it returns FALSE;

The x value returned by snaprounding "seems" to be as close as possible from 311.4 SELECT ST_X(ST_GeomFromText('POINT(0 311.4)')); -> display 311.4

But while SELECT ST_X(ST_GeomFromText('POINT(311.4 0)'))-311.4 display 0 SELECT ST_X(ST_SnapToGrid(ST_GeomFromText('POINT(311.4 0)'), 0.1))-311.4 display 5.6843418860808e-014

x-value is not the same after snap rounding.

My use case was : I want to check that all my coordinates are decimetric.

[pramsey]

Unfortunately this seems a little too deep to do anything about. We're seeing here the fact that decimal floats don't always have exact analogues in IEEE double precision (binary) storage. Here's a simple C program that demonstrates exactly the effect you're seeing.

static void test_rint()
{
    char str[] = "311.4";
    double d1, d2, d3;
    double g = 0.1;

    sscanf(str, "%lf", &d1);
    printf("d1: %g\n", d1);
    
    d2 = rint(d1 / g) * g;
    printf("d2: %g\n", d2);
    printf("d2-d1: %g\n", d2-d1);

    d3 = rint(d2 / g) * g;    
    printf("d3-d2: %g\n", d3-d2);
}

The results are

d1: 311.4
d2: 311.4
d2-d1: 5.68434e-14
d3-d2: 0

Basically the way scanf is forming 311.4 is slightly different from the version that the math operations of dividing and then remultiplying the number come up with. In order to "snap to a grid" all we're doing is scaling up numbers until the grid is represented by integral values, then rounding, so exactly the operation you're seeing in this example.

[michaudm]

Hi Paul,

I'm not a C programmer so I don't know if the problem can be solved in Geos/PostGIS, but in java,

rint(311.4/0.1)*0.1 -> 311.40000000000003

while

round(311.4/0.1)*0.1 -> 311.4

If I use rint, I obtain exactly the same result as in your example while with round, I get correct result.

Moreover, I noticed the following comment in the makePrecise method of JTS code :

if (modelType == FIXED) {
    return Math.round(val * scale) / scale;
    // return Math.rint(val * scale) / scale;
}

Not sure how to interpret these results as I'm not a IEEE double precision expert, but maybe there is something to do, especially if it makes JTS and Geos behaviour closer.

[pramsey]

Good thought, but it seems to not make the difference we want in C world.

static void test_rint()
{
    char str[] = "311.4";
    double d1, d2, d3;
    double g = 0.1;

    sscanf(str, "%lf", &d1);
    printf("d1: %g\n", d1);
    
    d2 = rint(d1 / g) * g;
    printf("d2: %g\n", d2);
    printf("d2-d1: %g\n", d2-d1);

    d3 = round(d1 / g) * g;    
    printf("d3-d1: %g\n", d3-d1);
}

d1: 311.4
d2: 311.4
d2-d1: 5.68434e-14
d3-d1: 5.68434e-14