I'm trying to calculate the departing angle of a geodesic using several different methods in two different coordinate spaces.
Currently I'm testing the following geodesic terminal points:
- Start (80°W, 40°N)
- End (40°E, 21°N)
The two different coordinate spaces are:
- Earth-ellipsoidal (WGS-84, $a = 6,378,137$, $f = 0.0033528106647475126$)
- Earth-spherical ($a = 6,378,137$)
The methods are:
- Plot the geodesic as a straight line in an orthographic projection (spherical only), and visually measure departing heading
- Infinitesimal departure distance from start, use
pyproj.Geod.inv_intermediateto generate start portion of geodesic path, apply
$$ \tan^{-1} \frac {\Delta y}{\Delta x} $$ That angle is clockwise from north and agrees very closely with method 1, producing 61.3° in the ellipsoidal and 61.5° in the spherical coordinate spaces respectively.
Green is method 1 and red is method 2.
The following methods agree exactly with each other but not the methods above:
pyproj.Geod.inv(), elliptical: 54.511°inv()spherical: 54.600°- Charles Karney's GeodSolve, elliptical: 54.511°
- GeodSolve, spherical: 54.600°
The GeodSolve method is described in C. F. F. Karney, Algorithms for geodesics.
Intuitively I would think that every result should be close, but the gap between the arctan method and the others is too large. I would prefer to use inv(), but my current assumptions are that
- GeodSolve implementation is correct
- Both
pyprojimplementations are "correct", but produce different quantities that have each have a different geometric definition - I've misunderstood what all of the theory is talking about when it describes departing angle
So what have I missed? How do I take the results from the second group of methods and match them up with the results of method 1?