[postgis-devel] ST_Project inaccuracy? Values out by many miles...

[Paul Ramsey]

Brent,

Had me going there for a sec, but Phew. It’s not me, it’s you.

Ok, here’s my thought experiment, and please forgive me for doing it in the northern hemisphere.

- I start at YVR and point my plane directly east, and start flying a great circle route, such that if I go on in a straight line, I’ll eventually hit YVR again. After 200 km am I (a) south (b) north (c) directly east of YVR?

Here’s one thing I know for sure, it’s not (c), because a line of latitude (except for the equator) is not a great circle. (It’s (a).)

Or, let’s do another one.

- You stand here at latitude 45. I’ll stand 200 miles to the east, on latitude 45. What bearing will you take to get to me, using the shortest possible path?

select degrees(ST_Azimuth('POINT(-126 45)'::geography, 'POINT(-125 45)'::geography));

Seeing it? Stand further away.

select degrees(ST_Azimuth('POINT(-126 45)'::geography, 'POINT(-25 45)'::geography));

The world is not flat and flat reasoning will lead you to the wrong places.

In your case, going directly east one way and then directly west the other, should lead to a result to the north of where you started, (a zig-zaggin path) as indeed it does.

P.

-- 
Paul Ramsey
http://cleverelephant.ca
http://postgis.net

On December 16, 2013 at 2:00:18 PM, Brent Wood (pcreso at pcreso.com) wrote:




Hi,

I have been looking at using ST_Project() to solve a problem I have (& in theory it is perfect!).

It seems to introduce some inaccuarcies, which can be quite significant in some use cases. The SQL below projects a new point 2km E, then back 2km W. The result should be in the same place we started, but is not, by a larger value than I'm comfortable with. The source point location comes from a real dataset.

select Postgis_full_version();

POSTGIS="2.0.1 r9979" GEOS="3.3.8-CAPI-1.7.8" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.2, released 2012/10/08" LIBXML="2.7.8" TOPOLOGY RASTER


select ST_AsText(point) as point,
       ST_AsText(ST_Project(point, 2000,radians(90))) as proj2k,
       ST_AsText(ST_Project(ST_Project(point, 2000,radians(90)), 2000,radians(270))) as proj2k0
from (select ST_AsText(ST_SetSRID(ST_MakePoint(171.0259489718, -45.2961933568),4326))::geography as point) as foo;

point   | POINT(171.0259489718   -45.2961933568)
proj2k  | POINT(171.051446314948 -45.2961905108324)
proj2k0 | POINT(171.025948973075 -45.296187664865)

The final lat & lon should be the same as the original, but both have changed, the latitude substantially more than the longitude - even though the point latitude should be unchanged (it only moved E & W, not N or S) 

It gets worse with larger distances, eg, at 2000km:

select ST_AsText(point) as point,
       ST_AsText(ST_Project(point, 2000000,radians(90))) as proj2000k,
       ST_AsText(ST_Project(ST_Project(point, 2000000,radians(90)), 2000000,radians(270))) as proj2000k0
from (select ST_AsText(ST_SetSRID(ST_MakePoint(171.0259489718, -45.2961933568),4326))::geography as point) as foo;

point      | POINT(171.0259489718 -45.2961933568)
proj2000k  | POINT(-164.264848338683 -42.5385273590618)
proj2000k0 | POINT(172.015779821441 -40.0225960863261)

it causes over 5 degree change in the latitude.

Is this a bug I should be submitting or have I missed something?


Thanks

  Brent Wood



_______________________________________________  
postgis-devel mailing list  
postgis-devel at lists.osgeo.org  
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-devel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-devel/attachments/20131216/d27d82f1/attachment.html>