Difference between Vincenty and great-circle distance calculations?

Question

Python's geopy package features two distance measurements techniques: Great Circle and Vincenty's formulae.

>>> from geopy.distance import great_circle
>>> from geopy.distance import vincenty
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> vincenty(p1, p2).meters
429.16765838976664
>>> great_circle(p3, p4).meters
428.4088367903001

What is the difference? Which distance measurement is preferred?

Answer 1

According to Wikipedia, Vincenty's formula is slower but more accurate:

Vincenty's formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of a spheroid, developed by Thaddeus Vincenty (1975a) They are based on the assumption that the figure of the Earth is an oblate spheroid, and hence are more accurate than methods such as great-circle distance which assume a spherical Earth.

The accuracy difference is ~0.17% in a 428 meters distance in Israel. I've made a quick-and-dirty speed test:

<class 'geopy.distance.vincenty'>       : Total 0:00:04.125913, (0:00:00.000041 per calculation)
<class 'geopy.distance.great_circle'>   : Total 0:00:02.467479, (0:00:00.000024 per calculation)

Code:

import datetime
from geopy.distance import great_circle
from geopy.distance import vincenty
p1 = (31.8300167,35.0662833)
p2 = (31.83,35.0708167)

NUM_TESTS = 100000
for strategy in vincenty, great_circle:
    before = datetime.datetime.now()
    for i in range(NUM_TESTS):
        d=strategy(p1, p2).meters
    after = datetime.datetime.now()
    duration = after-before
    print "%-40s: Total %s, (%s per calculation)" % (strategy, duration, duration/NUM_TESTS)

To conclude: Vincenty's formula is doubles the calculation time compared to great-circle, and its accuracy gain at the point tested is ~0.17%.

Since the calculation time is negligible, Vincenty's formula is preferred for every practical need.

Update: Following the insightful comments by whuber and cffk's and cffk's answer, I agree that the accuracy gain should be compared with the error, not the measurement. Hence, Vincenty's formula is a few orders of magnitude more accurate, not ~0.17%.

Answer 2

If you're using geopy, then the great_circle and vincenty distances are equally convenient to obtain. In this case, you should almost always use the one that gives you the more accurate result, i.e., vincenty. The two considerations (as you point out) are speed and accuracy.

Vincenty is two times slower. But probably in a real application the increased running time is negligible. Even if your application called for a million distance calculations, we are only talking about a difference in times of a couple of seconds.

For the points you use, the error in vincenty is 6 μm and the error in the great circle distance is 0.75 m. I would then say that vincenty is 120000 times more accurate (rather than 0.17% more accurate). For general points, the error in the great circle distance can be as much as 0.5%. So can you live with a 0.5% error in distances? For casual use (what's the distance from Cape Town to Cairo?), probably you can. However, many GIS applications have much stricter accuracy requirements. (0.5% is 5m over 1km. That really does make a difference.)

Nearly all serious mapping work is carried out on the reference ellipsoid and it therefore makes sense that distances should be measured on the ellipsoid too. Maybe you can get away with great-circle distances today. But for each new application, you will have to check whether this is still acceptable. Better is just to use the ellipsoidal distance from the start. You'll sleep better at night.

ADDENDUM (May 2017)

In reply to the answer given by @craig-hicks. The vincenty() method in geopy does have a potentially fatal flaw: it throws an error for nearly antipodal points. The documentation in the code suggests increasing the number of iterations. But this is not a general solution because the iterative method used by vincenty() is unstable for such points (each iteration takes you further from the correct solution).

Why do I characterize the problem as "potentially fatal"? Because any use of the distance function within another software library needs to be able to handle the exception. Handling it by returning a NaN or the great-circle distance may not be satisfactory, because the resulting distance function will not obey the triangle inequality which precludes its use, e.g., in vantage-point trees.

The situation isn't completely bleak. My python package geographiclib computes the geodesic distance accurately without any failures. The geopy pull request #144 changes the geopy's distance function to use geographiclib package if it's available. Unfortunately this pull request has been in limbo since Augest 2016.

ADDENDUM (May 2018)

geopy 1.13.0 now uses the geographiclib package for computing distances. Here's a sample call (based on the example in the original question):

>>> from geopy.distance import great_circle
>>> from geopy.distance import geodesic
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> geodesic(p1, p2).meters
429.1676644986777
>>> great_circle(p1, p2).meters
428.28877358686776

Answer 3

My apologies for posting a second answer here, but I taking the opportunity to respond to the request by @craig-hicks to provide accuracy and timing comparisons for various algorithms for computing the geodesic distance. This paraphrases a comment I make to my pull request #144 for geopy which allows the use of one of two implementations of my algorithm for geodesics to be used within geopy, one is a native python implementation, geodesic(geographiclib), and the other uses an implementation in C, geodesic(pyproj).

Here is some timing data. Times are in microsecs per call

method                          dist    dest
geopy great_circle              20.4    17.1
geopy vincenty                  40.3    30.4
geopy geodesic(pyproj)          37.1    31.1
geopy geodesic(geographiclib)  302.9   124.1

Here is the accuracy of the geodesic calculations based on my Geodesic Test Set. The errors are given in units of microns (1e-6 m)

method                        distance destination
geopy vincenty                 205.629  141.945
geopy geodesic(pyproj)           0.007    0.013
geopy geodesic(geographiclib)    0.011    0.010

I've include hannosche's pull request #194 which fixes a bad bug in the destination function. Without this fix, the error in the destination calculation for vincenty is 8.98 meters.

19.2% of the tests cases failed with vincenty.distance (iterations = 20). However the test set is skewed towards cases which would case this failure.

With random points on the WGS84 ellipsoid, the Vincenty algorithm is guaranteed to fail 16.6 out of 1000000 times (the correct solution is an unstable fixed point of the Vincenty method).

With the geopy implementation of Vincenty and iterations = 20, the failure rate is 82.8 per 1000000. With iterations = 200, the failure rate is 21.2 per 1000000.

Even though these rates are small, failures can be quite common. For example in a dataset of 1000 random points (think the worlds airports, perhaps), computing the full distance matrix would fail on average 16 times (with iterations = 20).

Answer 4

Whether using vincenty or haversine or the spherical law of cosines, there is wisdom in becoming aware of any potential issues with the code you are planning to use, things to watch out for and mitigate, and how one deals with vincenty vs haversine vs sloc issues will differ as one becomes aware of each one's lurking issues/edgecases, which may or may not be popularly known. The seasoned programmer knows this. Newbies may not. I hope to spare some of them frustration when a snippet from a forum does something unexpected, in certain cases. If one is seriously going to use some version of any of these, vincenty, haversine, sloc, then SE, SO, Reddit, Quora, etc, may have provided limited help in some initial coding of a solution, but that does not mean that their solution or accepted 'answer' is free of issues. If a project is important enough, it deserves an appropriate reasonable amount of research. Read the manual, read the docs, and if a code review of that code exists, read that. Copying and pasting a snippet or gist that was upvoted a hundred or more times does not mean its safety is comprehensive and assured.

The intriguing answer posted by cffk raises the point of being aware of lurking edgecases, in packaged solutions, that can produce exceptions or other difficulties. The specific claims made in that post are beyond my time budget to pursue at present, but I take away from it that there are indeed lurking issues in certain packages, including at least one vincenty implementation, regarding which at least one person has proposed to improve one way or another, in order to minimize or eliminate the risk of encountering those difficulties. I won't add further to that topic concerning vincenty (being far too ignorant of it), but will turn instead to haversine, at least partly on topic with the OP.

The popularly published haversine formula, whether in python or another language, because it is going to be most likely using the IEEE 754 floating point spec on most all intel and intel-like systems today, and ARM processors, powerPC, etc, it is going to also be susceptible to rare but real and repeatable exception errors very near or at 180 degree arc distance, antipodal points, due to floating point approximations and rounding. Some newbies may not yet have been bitten by this situation. Because this fp spec approximates and rounds, this doesn't mean that any code that calls on fp64 might cause exception errors, no. But some code, some formulas might have not so obvious edgecases where the approximations and roundings of IEEE 754 fp64 may cause a value to stray slightly out of the domain of a math method that is expected to flawlessly evaluate such a value. An example... sqrt(). If a negative value finds its way into a sqrt(), such as sqrt(-0.00000000000000000122739), there will be an exception error. In the haversine formula, the manner in which it progresses towards a solution, there are two sqrt() methods in the atan2(). The a that is calculated and then used in the sqrt(), can, at the antipodal points on the globe, slightly stray below 0.0 or above 1.0, very slightly because of fp64 approximations and rounding, rarely, but repeatably. Consistent reliable repeatability, in this context, makes this an exception risk, an edgecase to protect, to mitigate, rather than an isolated random fluke. Here's an example of a short python3 snippet of haversine, without the necessary protection:

import math as m

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

Very near or at antipodal points, a calculated in the first line of the formula may stray negative, rarely, but repeatably with those same lat lon coordinates. To protect/correct those rare occurrences, one can simply add, after the a calculation,as seen below:

import math as m

note = ''

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
if a < 0.0: a = 0.0 ; note = '*'
if a > 1.0: a = 1.0 ; note = '**'
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

# note = '*'  # a went below 0.0 and was normalized back to 0.0
# note = '**' # a went above 1.0 and was normalized back to max of 1.0

Of course I did not show the entire function here, but a short snippet as is so often posted. But this one shows the protection for the sqrt(), by testing the a, and normalizing it if necessary, also saving the need to put the whole thing in a try except. The note = '' up top is to prevent the bytecode stage from protesting that note is being used before being assigned a value, if it is returned with the result of the function.

With this simple change, of adding the two a tests, the sqrt() functions will be happy, and the code now has an additional note that can be returned to calling code, to alert that a result has been slightly normalized, and why. Some may care, some may not, but its there, preventing an exception error, that 'can' otherwise occur. A try except block may catch the exception, but not fix it, unless explicitly written to do so. It seems easier to code the correction line(s) immediately after the a calculation line. Thoroughly scrubbed input should then not require a try except block here at all.

Summary, if using haversine, coded explicitly rather than using a package or library, no matter your language of choice, it would be a good idea to test and to normalize a back into the needful range of 0.0 <= a <= 1.0 in order to protect the next line with its c calculations. But the majority of haversine code snippets do not show it, and do not mention the risk.

Experience: during thorough testing around the globe, in 0.001 degree increments, I've filled up a hard drive with lat lon combinations that caused an exception, a reliable consistent repeatable exception, during a month of also collaterally testing the reliability of the CPU cooling fan, and my patience. Yes, I've since deleted most of those logs, since their purpose was mostly to prove the point (if the pun is allowed). But I have some shorter logs of 'problem lat lon values', kept for testing purposes.

Accuracy: Will a and the entire haversine result lose some accuracy by normalizing it that small bit back into domain? Not much, maybe no more than the fp64 approximations and roundings were already introducing, that caused that slight drift out of domain. If you have found haversine acceptable over vincenty already -- simpler, faster, easier to customize, troubleshoot and maintain, then haversine may be a good solution for your project.

I've used haversine on an overhead projected skysphere for measuring angular distances between objects in the sky, as viewed from a position on earth, mapping azimuth and alt to skysphere lat lon equivalent-like coordinates, no elipsoid to consider at all, since the projected theoretical skysphere is a perfect sphere, when it comes to measuring angular distance look angles between two objects from a position on the earth's surface. It suits my needs perfectly. So, haversine still is very useful, and very accurate, in certain applications (well within my purposes) ... but if you do use it, whether on the earth for GIS or navigation, or in sky object observations and measurements, do protect it in the case of antipodal points or very near antipodal points, by testing a and nudging it back into its needful domain when needed.

The unprotected haversine is all over the internet, and I've only seen one old usenet post that showed some protection, I think from someone at JPL, and that may have been pre-1985, pre- IEEE 754 floating point spec. Two other pages mentioned possible issues near antipodal points, but did not describe those issues, or how one might mitigate them. So there is concern for the newbies (like me) who may not always understand good practice well enough to further research, and test edgecases, of some code they have copied and pasted into a project in trust. cffk's intriguing post was refreshing in that it was public with these types of issues, that are not often mentioned, rarely publicly coded for protection in snippets, and rarely discussed in this way, compared to the amount of unprotected and undiscussed versions that are posted.

As of 20190923, the wiki page for haversine formula does indeed mention the problem possible at antipodal points, due to floating point issues in computing devices ... encouraging ...

https://en.wikipedia.org/wiki/Haversine_formula

(because that wiki page does not, at this time, have an html anchor for the section to which I would directly link, therefore, after the page loads, do a search on that browser page for 'When using these formulae' and you will see the haversine's problem with antipodal points mentioned, more officially.)

And this other site also has a very brief mention of it:

https://www.movable-type.co.uk/scripts/latlong.html

If one does a find on that page for 'including protection against rounding errors', there is this...

If atan2 is not available, c could be calculated from 2 ⋅ asin( min(1, √a) ) (including protec­tion against rounding errors).

Now there is a rare instance where rounding errors are mentioned, and protection shown for the asin() version, yet not mentioned or shown for the atan2() version. But at least the risk of rounding errors are mentioned.

imho, any 24/7/365 application using haversine, needs this protection near the antipodal points as an important and simple detail.

I don't know which haversine packages do or do not include this protection, but if you are new to all this, and you are going to use the popularly published 'snippet' version(s), now you know it needs protection, and that protection is very simple to implement, that is, if you are not using vincenty, and not using a packaged haversine without easy access to modify the code of the package.

IOW, whether using vincenty or haversine or sloc, one ought to become aware of any issues with the code, things to watch out for and mitigate, and how one deals with vincenty vs haversine vs sloc issues will differ as one becomes aware of each one's lurking issues/edgecases, which may or may not be popularly known.

Answer 5

It appears that the geopy.distance package offers a function "distance()" which defaults to vincenty(). I would recommend using distance() on principle, as it is the package recommendation, in case is ever diverged from vincenty() in the future (unlikely as that is). Continue reading:

This documentation note is included in the source code for the vincenty() function you specified:

Note: This implementation of Vincenty distance fails to converge for some valid points. In some cases, a result can be obtained by increasing the number of iterations (iterations keyword argument, given in the class __init__, with a default of 20). It may be preferable to use :class:.great_circle, which is marginally less accurate, but always produces a result.

The source code with th above comment/note can be found at https://github.com/geopy/geopy/blob/master/geopy/distance.py Scroll down to the definition for vincenty()

Nevertheless, the default distance function used by that package when caliing distance() is the vincenty() function, which implies that failure to converge is not catastrophic, and a reasonable answer is returned - most importantly an exception is not generated.

Update: As noted by "cffk", the vincenty() function does explicitly throw a ValueError exception when the algorithm does not converge - although it is not documented in the function description. Hence, the documentation is buggy.