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

[Paul Ramsey]

Just to finish off with a visual flourish, here's what going 2000 km
east, then turning around and going west looks like.

http://www.gcmap.com/mapui?P=48.5N+123.4W-97.369416581944W+45.4347105567823N-122.131748813688W+42.6643751114942N&MS=wls&DU=mi

And the one sentence explanation "on the sphere, straight lines (what
you get by moving forward and not turning) are not lines of constant
bearing (except for meridians and the equator)"

P.


On Mon, Dec 16, 2013 at 2:15 PM, Paul Ramsey <pramsey at cleverelephant.ca> wrote:
> 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