Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this 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: – user1158559 Mar 10 '13 at 2:15
@user1158559, John Snyder's Map Projections: A Working Manual should have it: – 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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.