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

[realityexists]

This was broken again in r10670. Looks like the call to lwgeom_nudge_geodetic() was lost in that commit. With that change the test cases above no longer raise errors, but instead wrap around the dateline: -180 is returned instead of 180. It's the same underlying problem, though. Please add a unit test for this.

[realityexists]

Very simple patch attached. I'd really appreciate if this could get reviewed and committed, as it's stopping us from upgrading and we really want to upgrade to get the benefit of the #1828 fix.

I couldn't get the regression tests running, but have tried to add one anyway (see #2066).

[realityexists]

Now with a proper, working, regression test

[robe]

patch looks fine to me. pramsey -- your call?

[strk]

The regression test isnt' clear to me. What is it checking ? ST_AsText output is known not to show discrepancies... Why not taking ST_SnapToGrid out of that test ? And why is the NOTICE gone ?

[realityexists]

It checks that the points at 180 longitude remain at 180. Without this patch they become -180 instead. It also checks multi-geometries and a polygon with multiple rings, which have failed in the past (at r8990). ST_SnapToGrid is what causes the coordinates to be "just over" 180, so it's essential to the test.

The NOTICE is gone from another test because it no longer occurs. Looking at that test, I think it never should have occurred and was just allowed because it didn't seem to cause any problems, but maybe someone more familiar with it can comment.

[pramsey]

Sorry if I'm being obtuse, but 2.0 seems to include nudging behavior already. Your patch applies cleanly to trunk (2.1) and I've done so at r11109. If everything is OK, feel free to close. Thanks for the patch.

[realityexists]

Thank you! Tested r11113 and all seems to be well.

The problem only occurred on trunk. 2.0 would have been fine. I've changed the Milestone on this ticket to 2.1.0 in case that caused any confusion.

[realityexists]

Sorry, there is still a bug. My patch doesn't handle the case where one coordinate needs to be "nudged", but the other needs to be "coerced" (from #799). In that case it only did the nudging and returned a geography with one coordinate out of range! Second patch attached (with test).

[realityexists]

Paul, could you apply the second patch? I feel bad about my bug still being in there. :)

[robe]

committed at r11345 pramsey can change if he has issues.