23

Bug introduced in 8.0.1 or earlier and fixed in 10.0.0

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
cffk
  • 561
  • 2
  • 13
  • 1
    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.'s missing motivation Jun 06 '13 at 03:56
  • 2
    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 06 '13 at 11:30
  • Interesting... if I crank up the setting of the option "InverseMaxIterations" for GeoDistance[], I still get the incorrect results for the first two, but I get something much nearer to the "correct" value for the third example. If you want to debug what's happening, you might want to set the option "InverseVerbose" to True. – J. M.'s missing motivation Jun 06 '13 at 11:51
  • 1
    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 06 '13 at 12:54
  • There's an alternative setting for Method, "ExtendedNewtonRaphson". It fails on your first example, but returns results close to the "correct" ones for the other two. – J. M.'s missing motivation Jun 06 '13 at 14:13
  • 1
    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 06 '13 at 15:48
  • Might I propose posting an implementation of an alternative algorithm (e.g. the improved Vincenty method) as an answer to your question? :) – J. M.'s missing motivation Jun 06 '13 at 15:54
  • With "I reported" you mean you sent a bug report to wolfram support? Or did you report it somewhere on a forum? – Sjoerd C. de Vries Jun 06 '13 at 21:36
  • 1
    I submitted it Wolfram support on 2011-06-27 and got a reply (TS 22876) saying it had been forwarded to "the appropriate people in our development group". I did a follow-up a few days later with a similar reply. That was the last I heard about the matter (and I couldn't pursue it further since my trial license ran out). – cffk Jun 07 '13 at 12:59
  • Maybe you could write a Mathematica implementation of the methods in this paper... :) – J. M.'s missing motivation Jun 08 '13 at 15:21
  • @J. M. the link is broken. – shrx Jul 01 '13 at 19:49
  • 3
    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 http://arxiv.org/abs/1109.4448 – cffk Jul 01 '13 at 20:00

3 Answers3

12

From the information provided by @J. M., 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)?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
cffk
  • 561
  • 2
  • 13
  • 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.'s missing motivation Jun 09 '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 09 '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.'s missing motivation Jun 09 '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.'s missing motivation Jun 09 '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 09 '13 at 17:19
11

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>
cffk
  • 561
  • 2
  • 13
  • 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 03 '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 04 '15 at 00:08
  • I know a lot of work went into making sure it was correct, but what exactly, I don't know. – rcollyer May 04 '15 at 02:20
  • 2
    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.'s missing motivation May 04 '15 at 04: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 04 '15 at 20:32
  • Dear J. M., I have no access to Mathematica at present. However, Carlson's symmetric integrals are definitely the way to go for solving geodesics (and I have this implemented in C++ with MPFR). – cffk May 04 '15 at 20:42
  • The "next few weeks" took longer than I would like, but I've finally implemented it. – J. M.'s missing motivation Mar 22 '19 at 10:31
8

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

David Creech
  • 1,128
  • 10
  • 15
  • 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 https://books.google.com/books?id=Zk0wAAAAIAAJ&pg=PA333 – cffk Jul 31 '15 at 19:20
  • I received a response from Wolfram Tech Support that they have forwarded an incident report to the developers regarding this. – David Creech Jul 31 '15 at 21:12
  • 2
    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 03 '15 at 12:40
  • 1
    @cffk Case (3) is fixed in 10.4.0. –  Aug 26 '16 at 16:30