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

Mon Dec 16 17:14:38 PST 2013

Yep, a ticket is the best place. In an ideal world, a patch on the docbook, but just words will do in a pinch.

P.

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

On December 16, 2013 at 5:12:22 PM, Brent Wood (pcreso at pcreso.com) wrote:

Cheers Paul,

That's cool. I appreciate the reasoning - it might be usefil if the ST_Project() docs explained this... I'm happy to provide new words - is there a process whereby  I can submit them for approval?

Thanks,

  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 1:54 PM
Subject: Re: [postgis-devel] ST_Project inaccuracy? Values out by many miles...

Hi Brent,

The idea is that ST_Project and ST_Azimuth and ST_Distance work in concert. You can take the azimuth and distance between two points, and use st_project to go from one to another. While you don’t find the behaviour helps your use case, imagine something like “I want to go 10% of the way to that point there”. Find the azimuth and distance, and project 10% of the distance in the required direction. 

Even a path of non-constant bearing has a *starting* bearing, and that’s what gets plugged into ST_project.

Best,

Paul

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

On December 16, 2013 at 4:48:55 PM, Brent Wood (pcreso at pcreso.com) wrote:
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/29f20477/attachment.html>