<html><body><div style="color:#000; background-color:#fff; font-family:Courier New, courier, monaco, monospace, sans-serif;font-size:12pt"><div><span style="font-family: arial,helvetica,sans-serif;"><span>Thanks Paul</span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><span style="font-family: arial,helvetica,sans-serif;"><span><br></span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><span style="font-family: arial,helvetica,sans-serif;"><span>I wondered if it was something I was missing in the spherical geography arena...</span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><br><span style="font-family:
 arial,helvetica,sans-serif;"><span></span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><span style="font-family: arial,helvetica,sans-serif;"><span>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!!</span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><br><span style="font-family: arial,helvetica,sans-serif;"><span></span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;">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.</div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><span style="font-family: arial,helvetica,sans-serif;"><span><br></span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><span style="font-family: arial,helvetica,sans-serif;"><span>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 :-(  </span></span>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. <br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;">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.<br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color:
 transparent; font-style: normal;"><br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;">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. <br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;">Cheers (& keep up the good work!)<br></div><div style="color: rgb(0, 0, 0); font-size: 16px;
 font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><br></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;">  Brent<br><span style="font-family: arial,helvetica,sans-serif;"><span></span></span></div><div style="color: rgb(0, 0, 0); font-size: 16px; font-family: arial,helvetica,sans-serif; background-color: transparent; font-style: normal;"><span style="font-family: arial,helvetica,sans-serif;"><span><br></span></span></div><div><br></div>  <div style="font-family: Courier New, courier, monaco, monospace, sans-serif; font-size: 12pt;"> <div style="font-family: times new roman, new york, times, serif; font-size: 12pt;"> <div dir="ltr"> <hr size="1">  <font face="Arial" size="2"> <b><span style="font-weight:bold;">From:</span></b> Paul Ramsey <pramsey@cleverelephant.ca><br> <b><span style="font-weight:
 bold;">To:</span></b> Brent Wood <pcreso@pcreso.com>; PostGIS Development Discussion <postgis-devel@lists.osgeo.org> <br> <b><span style="font-weight: bold;">Sent:</span></b> Tuesday, December 17, 2013 11:15 AM<br> <b><span style="font-weight: bold;">Subject:</span></b> Re: [postgis-devel] ST_Project inaccuracy? Values out by many miles...<br> </font> </div> <div class="y_msg_container"><br><div id="yiv5173736304"><style>#yiv5173736304 body{font-family:Helvetica, Arial;font-size:13px;}</style><div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;">Brent,</div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;">Had me going there for a sec, but Phew. It’s not me, it’s you.</div><div
 id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;">Ok, here’s my thought experiment, and please forgive me for doing it in the northern hemisphere.</div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;">- 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?</div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont"
 style="font-family:Helvetica, Arial;font-size:13px;margin:0px;">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).)</div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;">Or, let’s do another one.</div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica, Arial;font-size:13px;margin:0px;">- 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?</div><div id="yiv5173736304bloop_customfont" style="font-family:Helvetica,
 Arial;font-size:13px;margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="margin:0px;">select degrees(ST_Azimuth('POINT(-126 45)'::geography, 'POINT(-125 45)'::geography));</div><div id="yiv5173736304bloop_customfont" style="margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="margin:0px;">Seeing it? Stand further away.</div><div id="yiv5173736304bloop_customfont" style="margin:0px;"><br clear="none"></div><div id="yiv5173736304bloop_customfont" style="margin:0px;">select degrees(ST_Azimuth('POINT(-126 45)'::geography, 'POINT(-25 45)'::geography));</div> <div class="yiv5173736304bloop_sign" id="yiv5173736304bloop_sign_1387231571592606976"><div><br clear="none"></div><div>The world is not flat and flat reasoning will lead you to the wrong places.</div><div><br clear="none"></div><div>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.</div><div><br clear="none"></div><div>P.</div><div><br clear="none"></div><span style="font-family:helvetica, arial;font-size:13px;"></span>-- <br clear="none">Paul Ramsey<br clear="none">http://cleverelephant.ca<div>http://postgis.net</div></div> <br clear="none"><div class="yiv5173736304yqt0868786921" id="yiv5173736304yqt72713"><div style="color:#A0A0A8;">On December 16, 2013 at 2:00:18 PM, Brent Wood (<a rel="nofollow" shape="rect" ymailto="mailto://pcreso@pcreso.com" target="_blank" href="mailto://pcreso@pcreso.com">pcreso@pcreso.com</a>) wrote:</div> <blockquote class="yiv5173736304clean_bq" type="cite"><span></span><div><div>



</div></div></blockquote></div></div><div class="yiv5173736304yqt0868786921" id="yiv5173736304yqt06751"><title></title><div><div style="color:#000;background-color:#fff;font-family:Courier New, courier, monaco, monospace, sans-serif;font-size:12pt;">
<div><span><br clear="none"></span></div>
<br clear="none">
<div style="font-family:Courier New, courier, monaco, monospace, sans-serif;font-size:12pt;">
<div style="font-family:times new roman, new york, times, serif;font-size:12pt;">
<div class="yiv5173736304y_msg_container"><br clear="none">
<div id="yiv5173736304">
<div>
<div style="color:#000;background-color:#fff;font-family:times new roman, new york, times, serif;font-size:12pt;">
<div>Hi,</div>
<div><br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:times new roman, new york, times, serif;background-color:transparent;font-style:normal;">
I have been looking at using ST_Project() to solve a problem I have
(& in theory it is perfect!).</div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:times new roman, new york, times, serif;background-color:transparent;font-style:normal;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:times new roman, new york, times, serif;background-color:transparent;font-style:normal;">
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.<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:times new roman, new york, times, serif;background-color:transparent;font-style:normal;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:times new roman, new york, times, serif;background-color:transparent;font-style:normal;">
<span style="font-family:Courier New, courier, monaco, monospace, sans-serif;">select
Postgis_full_version();<br clear="none">
<br clear="none">
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<br clear="none"></span><br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:times new roman, new york, times, serif;background-color:transparent;font-style:normal;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<span style="font-family:Courier New, courier, monaco, monospace, sans-serif;">select
ST_AsText(point) as point,<br clear="none">
       ST_AsText(ST_Project(point,
2000,radians(90))) as proj2k,</span></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:times new roman, new york, times, serif;background-color:transparent;font-style:normal;">
<span style="font-family:Courier New, courier, monaco, monospace, sans-serif;">  
    ST_AsText(ST_Project(ST_Project(point,
2000,radians(90)), 2000,radians(270))) as proj2k0<br clear="none">
from (select ST_AsText(ST_SetSRID(ST_MakePoint(171.0259489718,
-45.2961933568),4326))::geography as point) as foo;</span></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:Courier New, courier, monaco, monospace, sans-serif;background-color:transparent;font-style:normal;">
<span style="font-family:Courier New, courier, monaco, monospace, sans-serif;"><br clear="none">
</span></div>
<div style="color:rgb(0, 0, 0);font-size:16px;font-family:Courier New, courier, monaco, monospace, sans-serif;background-color:transparent;font-style:normal;">
<span style="font-family:Courier New, courier, monaco, monospace, sans-serif;">point  
| POINT(171.0259489718   -45.2961933568)<br clear="none">
proj2k  | POINT(171.051446314948 -45.2961905108324)<br clear="none">
proj2k0 | POINT(171.025948973075 -45.296187664865)<br clear="none">
<br clear="none"></span></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;">
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) </div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
It gets worse with larger distances, eg, at 2000km:</div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<span style="font-family:Courier New, courier, monaco, monospace, sans-serif;">select
ST_AsText(point) as point,<br clear="none">
       ST_AsText(ST_Project(point,
2000000,radians(90))) as proj2000k,<br clear="none">
      
ST_AsText(ST_Project(ST_Project(point, 2000000,radians(90)),
2000000,radians(270))) as proj2000k0<br clear="none">
from (select ST_AsText(ST_SetSRID(ST_MakePoint(171.0259489718,
-45.2961933568),4326))::geography as point) as foo;<br clear="none">
<br clear="none">
point      | POINT(171.0259489718
-45.2961933568)<br clear="none">
proj2000k  | POINT(-164.264848338683 -42.5385273590618)<br clear="none">
proj2000k0 | POINT(172.015779821441 -40.0225960863261)</span></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:Courier New, courier, monaco, monospace, sans-serif;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:Courier New, courier, monaco, monospace, sans-serif;">
<span style="font-family:times new roman, new york, times, serif;">it causes
over 5 degree change in the latitude.</span></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
Is this a bug I should be submitting or have I missed
something?</div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
Thanks</div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
  Brent Wood<br clear="none"></div>
<div style="color:rgb(0, 0, 0);font-size:16px;background-color:transparent;font-style:normal;font-family:times new roman, new york, times, serif;">
<br clear="none"></div>
</div>
</div>
</div>
<br clear="none">
<br clear="none"></div>
</div>
</div>
</div>


_______________________________________________
<br clear="none">postgis-devel mailing list
<br clear="none">postgis-devel@lists.osgeo.org
<br clear="none">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-devel</div></div></div><br><br></div> </div> </div>  </div></body></html>