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

Mon Dec 16 17:12:21 PST 2013

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/1ffdf216/attachment.html>