GeographicLib::Geodesic and GeographicLib::GeodesicLine provide accurate solutions to the direct and inverse geodesic problems. The Geod utility provides an interface to these classes.
References
The shortest path between two points on the ellipsoid at (lat1, lon1) and (lat2, lon2) is called the geodesic. Its length is s12 and the geodesic from point 1 to point 2 has azimuths azi1 and azi2 at the two end points. (The azimuth is the heading measured clockwise from north. azi2 is the "forward" azimuth, i.e., the heading that takes you beyond point 2 not back to point 1.)
Given lat1, lon1, azi1, and s12, we can determine lat2, lon2, azi2. This is the direct geodesic problem. (If s12 is sufficiently large that the geodesic wraps more than halfway around the earth, there will be a true geodesic between the points with a smaller s12.)
Given lat1, lon1, lat2, and lon2, we can determine azi1, azi2, s12. This is the inverse geodesic problem. Usually, the solution to the inverse problem is unique. In cases where there are multiple solutions (all with the same s12, of course), all the solutions can be easily generated once a particular solution is provided.
On a sphere, the geodesic is a great circle. On the ellipsoid, the problem is usually solved by transferring the geodesic to an "auxiliary sphere" where the latitude phi has been replaced by the reduced latitude beta where tan(beta) = (1 - f) tan(phi). On this sphere, the geodesic is again a great circle and the azimuth is the same as on the ellipsoid. However, the ellipsoidal distance is related to the great circle distance by an integral; and the ellipsoidal longitude is similarly related to the longitude on the auxiliary sphere. Williams (2002) gives the necessary integrals.
A popular method of doing geodesic calculations is due to Vincenty (1975). This entails expanding the integrals in powers of the ellipsoid flattening, and then solving both the direct and inverse problems by iteration. This approach has three drawbacks:
Here we reformulate the geodesic problem with three goals:
We derive higher-order series for the distance and longitude integrals using Maxima. In this release, we use an 8th-order expansion. The series is evaluated by a combination of Horner's method and Clenshaw summation. For the direct problem, we use a reversion of the series for the distance to compute the spherical arc length and thereby complete the solution without the need to iterate.
For the inverse problem, we rearrange the points so that the point 1 is the one closest to a pole and is in the southern hemisphere and the point 2 is to the north and east of point 1. We consider geodesics leaving point 1 with azimuth alpha1 in [0, pi] and we determine the longitude where these geodesics first intersect the latitude = lat2. The longitude is an increasing function of alpha1 and we use Newton's method to determine the value of alpha1 for which the longitude matches lon2. This then gives the desired geodesic.
A test set a geodesics is available at
This is about 25 MB (compressed). This consists of a set of geodesics for the WGS84 ellipsoid.
Each line of the test set gives 7 space delimited numbers
These are computed using as direct geodesic calculations with the given lat1, lon1, azi1, and s12. The distance s12 always corresponds to an arc length on the auxiliary sphere of pi or less, so the given geodesics give the shortest paths from point 1 to point 2. For simplicity and without loss of generality, lat1 is chosen in [0o, 90o], lon1 is taken to be zero, azi1 is chosen in [0o, 180o]. Furthermore, lat1 and azi1 are taken to be multiples of 10-12 deg and s12 is a multiple of 0.1 um in [0 m, 20003931.4586254 m]. This results lon2 in [0o, 180o] and azi2 in [0o, 180o].
The direct calculation uses an expansion of the geodesic equations accurate to f20 (approximately 1 part in 1050) and is computed with with Maxima's bfloats and fpprec set to 100 (so the errors in the data are probably 1/2 of the values quoted above).
The contents of the file are as follows:
(a total of 500000 entries). The values for s12 for the geodesics running between vertices are truncated to a multiple of 0.1 pm and this is used to determine point 2.
This data can be fed to the Geod utility as follows
gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,3,7 | ./Geod
gunzip -c GeodTest.dat.gz | cut -d' ' -f4,5,6,7 | sed "s/ \([^ ]*$\)/ -\1/" | ./Geod
gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,4,5 | ./Geod -i
Add, e.g., "-p 6", to the call to Geod to change the precision of the output. Adding "-f" causes Geod to print all 7 fields specifying the geodesic (in the same order as the test set).
Another much smaller set of test points has been culled from the literature on geodesics. These are all for the International Ellipsoid (a = 6378388 m, f = 1/297.0). The endpoints are taken from
However some of the results quoted for the inverse problem in these papers are inaccurate. Accurate data are given here. The end points are
No lat1 lat2 lon2 - lon1 00 37d19'54.95367" 26d07'42.83946" 41d28'35.50729" 01 35d16'11.24862" 67d22'14.77638" 137d47'28.31435" 02 1d00'00" -0d59'53.83076" 179d17'48.02997" 03 1d00'00" 1d01'15.18952" 179d46'17.84244" 04 41d41'45.88" -41d41'46.20" 179d59'59.44" 05 0d00'00" 0d00'00" 179d41'49.78063" 06 30d00'00" -30d00'00" 179d40'00" 07 60d00'00" -59d59'00" 179d50'00" 08 30d00'00" -29d50'00" 179d48'00" 09 30d00'00" -29d55'00" 179d48'00" 10 34d28'44.76421" -34d28'44.76421" 179d30'00" 11 0d00'00" 0d00'00" 179d23'38.18182" 12 34d28'44.764213" -34d28'44.764213" 179d30'00" 13 56d41'58.297496" -56d41'58.297496" 179d40'00" 14 29d45'00" 29d45'00.23848" 0d00'00.25626" 15 41d41'45.88" 41d41'46.2" 0d00'00.56" 16 46d00'00" 46d00'01" 0d00'01.816" 17 40d00'00" 40d00'02" 0d00'04.75" 18 38d00'00" 38d00'04.765" 0d00'05.554" 19 30d00'00" 37d53'32.46584" 116d19'16.68843" 20 30d19'54.95367" -30d11'50.15681" 179d58'17.84244" 21 0d39'49.12586" -0d45'14.13112" 179d58'17.84244" 22 0d00'54.95367" 0d00'42.83946" 179d28'17.84244" 23 40d00'00" -40d00'05.75932" 179d55'15.59578" 24 37d00'00" 28d15'36.69535" 2d37'39.52918" 25 38d30'45" -35d25'35" 179d45'00" 26 60d00'00" -60d00'00" 179d41'47"
and the geodesics are given by
No azi1 azi2 s12 00 95d27'59.63088905553491" 118d05'58.96160858886728" 4085966.7025902201825 01 15d44'23.74849770324814" 144d55'39.92147266777383" 8084823.8382961415712 02 88d59'59.99897053696689" 91d00'06.11835637659787" 19959999.9998034962728 03 4d59'59.9999565312182" 174d59'59.88480004914149" 19780006.5587880182731 04 179d58'49.16247860972594" 0d01'10.83761892414861" 20004566.7228054132931 05 30d00'00.00002084815889" 149d59'59.99997915184111" 19996147.416826781925 06 39d24'51.80601183179472" 140d35'08.19398816820528" 19994364.6068583984182 07 29d11'51.07006518417487" 150d49'06.86792886226858" 20000433.9629039632049 08 16d02'28.33895348073478" 163d59'10.33689436714454" 19983420.1535833515076 09 18d38'12.55689701007199" 161d22'45.43724069005652" 19992241.7634404403113 10 89d59'59.78608904242135" 90d00'00.21391095757865" 19981603.2781440234735 11 89d59'51.5863776476884" 90d00'08.4136223523116" 19970827.8695289752144 12 90d00'00.00000000479721" 90d00'00.00000000479721" 19981603.2781440234735 13 90d00'00.0000006754891" 90d00'00.0000006754891" 19994364.6068583984183 14 43d09'29.0634310812721" 43d09'29.19059165112047" 10.0665488514951 15 52d40'39.39067110974434" 52d40'39.76317180931766" 16.2839750636094 16 51d41'11.54166708300919" 51d41'12.8479912183896" 49.8037586266135 17 61d18'00.21784992705284" 61d18'03.27110871411274" 128.4581417556041 18 42d41'10.35602545797828" 42d41'13.77545984493403" 199.8717376568129 19 45d00'00.00000438165649" 129d08'12.32600911217221" 10002499.999860115911 20 2d23'52.10812966032674" 177d36'19.67010864697937" 19989590.5480170316779 21 177d39'39.01010365081902" 2d20'21.15303599834387" 19994529.4454322340309 22 54d08'27.73161935117439" 125d51'32.27232710890937" 19977290.7711390609949 23 170d15'10.8812277902834" 9d44'49.94565915545135" 20003827.8511392345191 24 164d59'59.99997936397788" 166d25'16.2593819478948" 1000000.0001515200393 25 3d22'19.5694402667547" 176d45'41.43721792855263" 19661438.0251956080829 26 90d00'00.00891264563132" 90d00'00.00891264563132" 19996104.3689008382234
We regard all the test data as "exact". For each member of the test set, 7 error distances are computed.
To guard against a canceling error, the calculations in (2), (4), (6), and (7) use a long double version of the direct geodesic calculation (which is about 2000 times more accurate than the double version). The maximum of the resulting errors in less than 12 nm for the double version and 6 pm for the long double version.
None of the checks listed above directly compares the computed azimuths with the exact values. The problem is that tiny changes in the positions of the points may sometimes lead to large changes in azimuths. Obvious examples of this are when the two points are close to opposite poles or when they are very close to each other.
A less obvious case is when the latitudes of the two points are equal and opposite and the azimuths are close to +/- 90o. For example, consider to points on the WGS84 ellipsoid
The two endpoints 2a and 2b are about 0.4nm apart and consequently the two geodesics (from 1 to 2a and from 1 to 2b) are very nearly the same length. However the azimuths at point 1 are
and the two geodesics are separated by about 160 m at their midpoints.
So instead of directly quoting the errors in the azimuths, we must use the back-to-front formulation: The azimuths returned are the correct azimuths for some geodesic with end-points within 12 nm of the given end-points.
When either point is at a pole, the azimuth is defined by keeping the longitude fixed and writing lat = 90 - eps or -90 + eps and taking the limit eps -> 0 from above.
The following prescription allows you to ascertain when there are multiple solutions to the inverse problem and to generate the complete set.
A rather large fraction of the code in GeographicLib::Geodesic deals merely with the spherical trigonometry required to do great circle calculations. When high accuracy is required, care has to be taken to use stable methods for computing angles. Typically, we represent angles internally with the pair [sin(theta), cos(theta)]. This expands the number of representable angles by about 4; in particular, it allows angles very close to all the cardinal points to be resolved. We use expressions for sin(theta) and cos(theta) which avoid large cancellations. Having a comprehensive test set as described in Test data for geodesics helps to identify problems. This representation allows input arguments which are multiples of 90o to be represented exactly.
The code for the inverse treats equatorial and meridional geodesics specially. All other case are handled by the same code (implementing Newton's method) with tests to select appropriate starting guess in difference regimes. Specifying a zero inverse flattening in Geodesic constructor is equivalent to zero flattening (i.e., a sphere). The method accommodates this case without any special coding.
In order to limits the enumeration of cases, the inverse problem is brought to a canonical form by interchanging the points and changing the signs of the coordinates. The effect of these changes are undone at the end by applying sign changes to the final calls to atan2 to give the azimuths.
The evaluation of the derivative for Newton's method for the inverse problem is a straightforward application of the chain rule for derivatives. In evaluating the derivative of the longitude function, the derivative with respect to cos2 alpha is carried out by differentiating the series. However the derivative with respect to sigma (the limit of the integral) just gives the integrand.
Geographic coordinates which are close to zero are quantized (by adding and subtracting 1/16). This allows certain cases of underflow to be avoided. This limits the accuracy to which coordinates can be specified to 0.7 pm. This should not be a problem since the overall accuracy of the algorithms is only 12 nm.
GeodesicLine::Position provides a way of doing a series of direct calculations for a single geodesic. This is about 2.5 faster than calling Geodesic::Direct. On a 2.6 GHz Intel machine (g++, version 4.3.0, -O3), calls to Geodesic::Direct take 1.1 us, calls to GeodesicLine::Position take 0.44 us, and calls to Geodesic::Inverse take 3.4 us.
Future work:
1.7.5