# How accurate is GeoDistance?

Two years ago I tried an evaluation copy of Mathematica. I reported various inaccuracies in the GeoDistance routine to which I never received any fix. Could someone check whether GeoDistance is any more accurate now? Thanks.

Here are the problems I reported on June 20, 2011:

                                     Mathematica   Correct     Error
GeoDistance[{30, 0}, {-30, 180}]  -> 19928486.7   20003931.5  -75444
GeoDistance[{0.4, 0}, {0.4, 180}] -> 19915473.0   19915472.0       1
GeoDistance[{0, 0}, {0.5, 179.5}] -> 19936426.5   19936288.6     138

-
For completeness: which reference frame are you using for the "correct" versions? By default, Mathematica uses "ITRF00" as the reference frame, and GeoDistance[] uses Vincenty's method (Method -> "Vincenty75") for distances. – J. M. Jun 6 '13 at 3:56
The "correct" results were calculated for the GRS80 ellipsoid (a = 6378137, f = 1/298.25722210088271), i.e., the reference ellipsoid used by ITRF00. Using WGS84 gives the same results (when rounded to the number of digits given here). – cffk Jun 6 '13 at 11:30
The "Vincenty75" method has known limitations; so if GeoDistance is still using this method, the results will be suspect for nearly antipodal points. Incidentally, Vincenty himself recognized the problems with his method and supplied a fix that went some way to curing the problems. – cffk Jun 6 '13 at 12:54
Thanks, this is good to know. However, it does make you wish that Mathematica would get rid of its three slightly flaky methods for GeoDistance and instead use one that just worked! – cffk Jun 6 '13 at 15:48
I hope this is only a temporary situation as we paid $3k so that anyone could access this paper. If you continue to have problems, you can download the preprint version of the paper from arxiv.org/abs/1109.4448 – cffk Jul 1 '13 at 20:00 ## 3 Answers From the information provided by @0x4A4D, the answer to my question is that GeoDistance has the same problems as it did 2 years ago. I implemented (in C++) some improved algorithms for geodesics in 2009 which are described in this paper. The benefits relative to GeoDistance are: • accurate to round-off for$|f| < 1/50$; • always converges; • computes ellipsoidal areas (and some other properties). Subsequently I've ported the algorithms to C, Fortran, Python, MATLAB, and JavaScript (see http://geographiclib.sf.net). The port to Java is in progress (it'll be ready in less than a month). I'm not in a position to write a Mathematica version because I no longer have access to Mathematica. However, if someone else wants to take the lead on this, I would be happy to assist. Alternatively, perhaps it's easier to interface Mathematica to the C or Java version (Java would be better since it presents a class interface)? - It hopefully shouldn't be too hard. OTOH, since Mathematica has explicit support for elliptic integrals, one might be able to use those instead of the truncated Fourier approximations in your paper. – J. M. Jun 9 '13 at 13:53 The formulation in terms of elliptic integrals is given here and implemented in the C++ library (no other languages at this point). For WGS84 with C++, the elliptic integral method is roughly 2 times less accurate (more accumulated round-off error) and 2 times slower. On the other hand, elliptic integrals are good for any flattening (I've checked$1/100 < b/a < 100$). – cffk Jun 9 '13 at 14:20 "roughly 2 times less accurate" - depends on the implementation, I'd say. I note that you're using the imaginary modulus formulation (the Earth's oblate); I'm sure you're aware that those expressions can be reformulated so that the moduli are real and within$[0,1)$, which should make for more stable evaluation. OTOH, the formulae also involve$D(\phi,k)\$ and the elliptic integral of the third kind, and those are indeed very difficult to deal with numerically, so I'd guess that's why you're not getting full accuracy from the elliptic integral formulation. – J. M. Jun 9 '13 at 14:27
If push comes to shove, it might be more expedient to reformulate in terms of Carlson's symmetric integrals; at least, in that formulation, computing the two special functions I mentioned earlier is less bothersome. – J. M. Jun 9 '13 at 14:31
I use Carlson's symmetric integrals; this is important because the expressions in terms of Legendre's integrals involve big cancellations. The use of an imaginary modulus is a non-issue once the integrals have been converted to Carlson's forms. I don't find the increased error too surprising. With the series solution, round-off can be very well controlled. It is difficult to maintain this accuracy with the more general formulation (however, I admit that I didn't try really hard). In any case, the accuracy is still very good -- about 25 nanometers for WGS84. – cffk Jun 9 '13 at 17:19

I've completed the implementation of my geodesic algorithms in Java. Documentation (including download information) is available at

http://geographiclib.sf.net/html/java

I gather that it's easy to call Java code from Mathematica. If someone would be kind enough to post instructions on how to do this, I would be happy to include this in the documentation.

This implementation uses series expanded to 6th order in the flattening. This gives results accurate to round-off (using double precision) for |f| < 1/50. There are two possible extensions which may be of interest to Mathematica users:

1. The solution in terms of elliptic integrals (to deal with highly eccentric ellipsoids of revolution); this is currently implemented in the C++ and Maxima versions.

2. The use of arbitrary precision arbitrary precision arithmetic; this is currently implemented in the Maxima version. (Since version 1.37, released on 2014-08-08, this capability is also available in the C++ library.)

ADDENDUM: The Java library is now available from Maven Central (since 2015-04). To use it include the following dependency in your pom.xml

<dependency>
<groupId>net.sf.geographiclib</groupId>
<artifactId>GeographicLib-Java</artifactId>
<version>1.42</version>
</dependency>

-
I just checked 10.1, and the results are accurate. If you have a copy of it, would you check it to see if it meets your needs? – rcollyer May 3 '15 at 13:35
This is good to hear. Unfortunately, I don't have Mathematica, so I can't check for myself. From the documentation it looks like the old Methods options, Vincenty75, Robbin61, ExtendedNewtonRaphson, have disappeared. Presumably the function is now using some reliable method? – cffk May 4 '15 at 0:08
I know a lot of work went into making sure it was correct, but what exactly, I don't know. – rcollyer May 4 '15 at 2:20
Dear cffk, I have some plans to implement the numerical evaluation of the Carlson integrals in Mathematica in the next few weeks; this might be useful to you for your accurate implementation of the geodesic distances. If you're interested, drop me a line. – J. M. May 4 '15 at 4:19
I believe that Mathematica 10 is using elliptic integrals to solve for the geodesic. However, a naive formulation gives errors for geodesics which graze a pole. To check if this is a problem, try GeoDistance and GeoDirection for (1) {40,0} to {40,179.9999999999} and (2) {40,0} to {89.9999999999,97}, (3) {89.9999999999,97} to {40,0}. The right answers are: (1) 11144873.397924414760 m, 0.00000000007766797 deg, (2) 5572436.698963568587 m, 0.00000000012982366 deg, and (3) 5572436.698963568587 m, -82.99999999991642394 deg. – cffk May 4 '15 at 20:32

This appears to be fixed. Mathematica v10 matches your "Correct" values for all three examples.

-
David, thank you for following up on this. Indeed, v 10.2 matches the "Correct" values as well. – MarcoB Jul 31 '15 at 17:36
Could you also check that GeoDistance and GeoDirection for (1) {40,0} to {40,179.9999999999} and (2) {40,0} to {89.9999999999,97}, (3) {89.9999999999,97} to {40,0}? The right answers are: (1) 11144873.397924414760 m, 0.00000000007766797 deg, (2) 5572436.698963568587 m, 0.00000000012982366 deg, and (3) 5572436.698963568587 m, -82.99999999991642394 deg. Thanks. – cffk Jul 31 '15 at 17:43
They are slightly different. (1) gives 11144873.39798774 m and 0.00000000007791086747713033 deg, (2) gives 5572436.698995231 m and 0.0000000001298175198786149 deg, and (3) gives 5572436.698995231 m and -82.9990873583041 deg. Max relative error on distance is 5E-10 %, on direction it is 0.3 % for example (1). – David Creech Jul 31 '15 at 18:14
Thanks for these results. A good way of assessing them is to do the direct geodesic calculation starting at point 1 with the computed distance and azimuth and seeing how far you are from point 2. The results are (1) 63 um, (2) 32 um, and (3) 77.9 m. The last result is disappointingly large. The possible problems in the algorithm are (a) poor numerical hygiene when doing spherical trigonometry and (b) a failure to use Cayley's transformation of the elliptic integral for the longitude see books.google.com/books?id=Zk0wAAAAIAAJ&pg=PA333 – cffk Jul 31 '15 at 19:20
For comparison, the errors using GeographicLib (the C++ version using double precision) are about 4 nm (nanometers). There's no reason why Mathematica shouldn't be able to match this. – cffk Aug 3 '15 at 12:40