Tell me more ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 2 down vote favorite
1

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 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 at 1:43
I opened an issue in the project: trac.osgeo.org/proj/ticket/211 – user1158559 Mar 10 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 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 at 23:51

Know someone who can answer? Share a link to this question via email, Google+, Twitter, or Facebook.

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.