Faced with a requirement to compute many (hundreds of) millions of geodesic distances in managed (.NET Framework) code, the Applications Prototype Lab developed a C# implementation of the non-iterative Andoyer-Lambert-Thomas method. We present the code and the results of numerical experiments comparing the accuracy of this method with the well-known iterative method of Vincenty.
Vincenty and Andoyer-Lambert-Thomas Methods
Vincenty presented his method for solving direct and inverse geodesic problems in 1975 (see here). His method has been widely used in geodesy because it is accurate to within 1 mm on the Earth spheroid. The method is iterative and requires special handling of (near-)antipodal points where standard iterations may not converge.
The series expansion for the geodesic length in terms of the flattening f was published (with errors) by Forsyth in 1895. In the mid-20th century Andoyer and Lambert published expansions to the first order in f. In 1970 Thomas published a detailed report with direct and inverse solutions to the second order in f (see here).
Esri Projection Engine C API provides an implementation of the iterative Vincenty method (see here).
Implementation and Numerical Experiments
We implemented the non-iterative Andoyer-Lambert-Thomas method in C# (the code is available here). In our implementation the geodesic distance is computed with relative accuracy f².
We generated 10,000 geographic points distributed uniformly between the latitudes 25°N and 50°N and the longitudes 65°W and 125°W (this area covers the 48 coterminous states of the USA). About 50 millions pair-wise geodesic distances were computed using both methods and the following statistics were generated for the non-iterative method:
Generated 10000 Fibonacci Points
Pairwise Distances: 49995000
Max Absolute Error: 0.037175 m
RMS Absolute Error: 0.011609 m
Avg Absolute Error: 0.008524 m
Max Relative Error: 0.025049 ppm
RMS Relative Error: 0.006656 ppm
Avg Relative Error: 0.004755 ppm
Speed of the non-iterative Andoyer-Lambert-Thomas method was more than 7 times faster than the iterative Vincenty method with an average absolute distance error of less than 1 cm for the 48 coterminous states.
- Uniformly distributed geographic points were generated using Fibonacci lattice method (see http://arxiv.org/abs/0912.4540 ).
- Computation of all pair-wise geodesic distances is an “embarrassingly parallel” problem which was implemented using Microsoft Task Parallel Library available in .NET Framework 4 (see http://msdn.microsoft.com/en-us/library/dd460693(VS.100).aspx )
- A differential geometry description of Earth geodesics can be found at MathWorld http://mathworld.wolfram.com/OblateSpheroidGeodesic.html ; for a good introduction see Wikipedia http://en.wikipedia.org/wiki/Geographic_distance .
Contributed by Lenny K.