5

All coordinates expressed in this question are (longitude, latitude), the default for proj4. London is at (-0.1460838, 51.5170986), and San Francisco is at (-122.419302, 37.775549). According to the proj4's geod program, the geographic bearing of San Francisco from London is -43.204268653666211719, and the great circle is 8624363.307m long.

To find out where San Francisco is on an Azimuthal Equidistant projection centred in London, I run the following:

cs2cs -I +proj=aeqd +lat_0=51.5170986 +lon_0=-0.1460838 +units=m +to +proj=lonlat
-122.419302 37.775549

And it gives me:

-5909365.29 6298188.68

Next, I transform these coordinates back to geographic coordinates:

cs2cs +proj=aeqd +lat_0=51.5170986 +lon_0=-0.1460838 +units=m +to +proj=lonlat
-5909365.29 6298188.68

And it gives me:

57d53'40.285"W 37d48'10.166"N 0.000

Which is a totally different location. Have I done something wrong to get this result, or is it a bug?

closed as off-topic by mkennedy, aldo_tapia, whyzar, Andre Silva, Midavalo Dec 14 '17 at 2:18

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This problem cannot or can no longer be reproduced. Changes to the system or to the asker's circumstances may have rendered the question obsolete, or the question does not include a procedure to enable potential answerers to reproduce the same symptoms. Such questions are off-topic as they are unlikely to help future readers, but editing them to include more details can lead to re-opening." – mkennedy, aldo_tapia, whyzar, Andre Silva, Midavalo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • I can confirm these results. It seems like you're doing everything right. – L_Holcombe Mar 9 '13 at 20:13
  • Thanks for confirming. I Narrowed it down to the inverse projection (elliptical) in <pre>src/PJ_aeqd.c</pre> - It uses the asin function to calculate the longitude, and its range is between <pre>-pi/2</pre> and <pre>pi/2</pre>. Anyone know where I can find a derivation for the projection between the Azimuthal Equidistant and Elliptical Geodetic coordinate systems? – user1158559 Mar 10 '13 at 1:43
  • I opened an issue in the project: trac.osgeo.org/proj/ticket/211 – user1158559 Mar 10 '13 at 2:15
  • 1
    @user1158559, John Snyder's Map Projections: A Working Manual should have it: pubs.er.usgs.gov/publication/pp1395 – mkennedy Mar 10 '13 at 6:33
  • Cheers. Turns out the implementation correctly follows John Snyder's book. John Snyder mentions in an example that a quadrant correction might be necessary, but he doesn't state how exactly to do that while defining the inverse projection equations. To be fair, he states that the equations he gives are only accurate enough up to a range of 800km, and we're talking 10,000km until the bug is encountered. To fix this bug, I'll need a derivation of those equations, starting with the forward equations. – user1158559 Mar 11 '13 at 23:51