ST_SnapToGrid returns a value out of range

[realityexists]

PostGIS 2.0.0 as of 2011-11-08

I've found some odd behaviour in ST_SnapToGrid, which looks like a bug. When I run

SELECT ST_SnapToGrid(ST_Transform(ST_GeomFromText('POINT(180 50)', 4269), 4326), 0.00001)::geography

an error is returned:

ERROR: Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

Calling ST_AsText() on the same geometry returns "POINT(180 50)", as expected. Interestingly, if I either decrease or even increase the precision the error does not occur. Eg. this works fine:

SELECT ST_SnapToGrid(ST_Transform(ST_GeomFromText('POINT(180 50)', 4269), 4326), 0.000001)::geography

It's only the value 0.0001 that seems to cause a problem.

[pramsey]

Actual, 0.00004 and 0.00007 also cause problems... but not the values in between. So a bona fide interesting case.

[realityexists]

Happens with latitude as well (2.0.0 r8577, Windows).

SELECT ST_SnapToGrid(ST_GeomFromText('POINT(-141 90)', 4326), 0.00001)::geography

=> ERROR: Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

Is there any workaround other than round-tripping everything though ST_AsText()? This is causing problems in my application that are difficult to track down, but I still want to use ST_SnapToGrid.

[pramsey]

OK, so this is tricky. Ignore the text for the moment. SnapToGrid? builds a virtual grid out from the origin and then snaps your coordinates to it. The net result is things move, for example {{{postgis20=# SELECT sT_astext(ST_SnapToGrid(ST_Transform(ST_GeomFromText('POINT(179.8 50)', 4269), 4326), 0.00007)::geography);

st_astext

---

POINT(179.79997 50.00002)

}}} Now, when you print out doubles, you only see a small portion of their total precision. The fact is that snapping even to 0.00001 is perturbing the coordinates slightly. And sometimes slightly positively. Now you might say "0.00001 fits evening into 180" and you'd be right, but that's in base 10, and the computer is working in base 2. So the grid formed by 0.00001 might actually snap 180 to something slightly larger than 180. This is legal and correct.

The question is what to do about it. One possibility would be to provide a function to wrap the coordinates to the geodetic bounds. This would have the effect of translating your perturbed points from 180plusabit to -180plusabit. But once you got into geographics it would all be the same, effectively.

I don't think we want the cast to do the wrap automatically (at least, people have told me that would be bad). Although I'm increasingly unsure of that.

[realityexists]

OK, I see - thanks for the detailed explanation! It totally makes sense from an implementation point of view, but totally doesn't from the user point of view. Agreed - it's tricky. The entire point of me calling ST_SnapToGrid is to "round" the coordinates - I would never expect it to increase the precision. Having the cast to geography automatically wrap might be a nice thing in general, but it wouldn't solve this particular problem*. I can't see any user expecting 'POINT(180 50)' to snap to 'POINT(-179.9999999 50)' or even to 'POINT(-180 50)' (for any snap size value). The only result a user would expect here is 'POINT(180 50)'. So I don't think that 180plusabit is "correct", it just isn't "incorrect enough" to matter in most cases.

Maybe ST_SnapToGrid() can err on the side of 180minusabit here somehow? Or maybe the cast to geography (and any other functions that check the range) can treat 180plusabit as equal to 180? For the latter to work the epsilon used would have to be greater than any imprecision that any snap size could ever create, so the former sounds safer to me.

*) Yes, for larger snap size values that don't fit evenly into 180 wrapping becomes a genuine user issue, not just an implementation issue. Eg.

ST_SnapToGrid(ST_GeomFromText('POINT(180 50)', 4326), 7)

returns 'POINT(182 50)' Should that fail or wrap when cast to geography? I'd say it should wrap, but that's probably outside the scope of this ticket.

[pramsey]

OK, going with a user-friendly approach, and when casting to geography am pushing points that are oh-so-close-but-not-quite at the boundaries over onto the boundaries. At 1e-10 hopefully this won't cause anyone any harm but will remove this minor problem. Particularly for snapping using tolerances that one would think would divide evenly into the bounds, this seems to be working admirable.

Trunk at r8980

[realityexists]

Thanks, I think that's a very good approach. Just tested it out in r8990.

This still fails:

SELECT ST_SnapToGrid(ST_GeomFromText('POINT(180 90)', 4326), 0.00001)::geography

This also fails:

WITH data
AS
(
	SELECT ST_GeomFromText('POLYGON((140 50,150 50,180 50,140 50),(140 60, 150 60,180 60,140 60))', 4326) AS one
	, ST_GeomFromText('GEOMETRYCOLLECTION(POINT(180 50),POINT(180 60))', 4326) AS two
)
SELECT ST_AsText(ST_SnapToGrid(one, 0.00001)::geography), ST_AsText(ST_SnapToGrid(two, 0.00001)::geography)
FROM data

[pramsey]

My code didn't take into account the possibility that a point could be out of bounds on two axes at once, which your brutal example exercised... fixed(?) at r8992

[realityexists]

Thanks, all good at r8992