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

Mon Dec 16 16:48:54 PST 2013

Thanks Paul

I wondered if it was something I was missing in the spherical geography arena...

I figured there would be some precision based error - so you would never quite get back to the same point - so ran the query to find this - then got slightly larger differences than I'd expected - so it was a very useful learning exercise!!

I'm still a bit confused - more about why this was done this way -  I don't think you can head 90 (due E) at a lat of 45 & be traversing a great circle (, especially without changing your bearing as you go) So the value of radians(90) passed to ST_Project() is not actually the direction to the new point. Which seems misleading.

The ST_Project docs don't mention great circles, or shortest distance, so I figured E (90) or W (270) just followed the constant line of longitude at that latitude, rather than disappearing on a great circle :-(  Also (with a few exceptions) a great circle is a line of constantly changing bearing - but the ST_Project() parameters specify a fixed azimuth value. 

I can see why the examples you gave return a non-90 degree heading - we are asking for the shortest path between two points - but the ST_Project() case is quite different - I'm traversing a course of 90 for a specified distance - which is not a great circle course.

That said, in my use case, the max distance is about 15nm, which gives a difference of < 1 degree at 45S for an E/W azimuth (which should be the largest difference)which is not enough to be an issue in this use case.

Cheers (& keep up the good work!)

  Brent

________________________________
 From: Paul Ramsey <pramsey at cleverelephant.ca>
To: Brent Wood <pcreso at pcreso.com>; PostGIS Development Discussion <postgis-devel at lists.osgeo.org> 
Sent: Tuesday, December 17, 2013 11:15 AM
Subject: Re: [postgis-devel] ST_Project inaccuracy? Values out by many miles...

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/33dda199/attachment.html>