To first approximation, Earth is a sphere. But it bulges at the equator, and to second approximation, Earth is an oblate spheroid. Earth is not exactly an oblate spheroid either, but the error in the oblate spheroid model is about 100x smaller than the error in the spherical model.

Finding the distance between two points on a sphere is fairly simple. Here’s a calculator to compute the distance, and here’s a derivation of the formula used in the calculator.

Finding the distance between two points on an ellipsoid is much more complicated. (A spheroid is a kind of ellipsoid.) Wikipedia gives a description of **Vincenty’s algorithm** for finding the distance between two points on Earth using an oblate spheroid model (specifically **WGS-84**). I’ll include a **Python** implementation below.

## Comparison with spherical distance

How much difference does it make when you calculate difference on an oblate spheroid rather than a sphere? To address that question I looked at the coordinates of several cities around the world using the `CityData`

function in Mathematica. Latitude is in degrees north of the equator and longitude is in degrees east of the prime meridian.

|--------------+--------+---------| | City | Lat | Long | |--------------+--------+---------| | Houston | 29.78 | -95.39 | | Caracas | 10.54 | -66.93 | | London | 51.50 | -0.12 | | Tokyo | 35.67 | 139.77 | | Delhi | 28.67 | 77.21 | | Honolulu | 21.31 | -157.83 | | Sao Paulo | -23.53 | -46.63 | | New York | 40.66 | -73.94 | | Los Angeles | 34.02 | -118.41 | | Cape Town | -33.93 | 18.46 | | Sydney | -33.87 | 151.21 | | Tromsø | 69.66 | 18.94 | | Singapore | 1.30 | 103.85 | |--------------+--------+---------|

Here are the error extremes.

The spherical model underestimates the distance from London to Tokyo by 12.88 km, and it overestimates the distance from London to Cape Town by 45.40 km.

The relative error is most negative for London to New York (-0.157%) and most positive for Tokyo to Sidney (0.545%).

**Update**: The comparison above used the same radius for both the spherical and ellipsoidal models. As suggested in the comments, a better comparison would use the *mean* radius (average of equatorial and polar radii) in the spherical model rather than the equatorial radius.

With that change the most negative absolute error is between Tokyo and New York at -30,038 m. The most negative relative error is between London and New York at -0.324%.

The largest positive error, absolute and relative, is between Tokyo and Sydney. The absolute error is 29,289 m and the relative error is 0.376%.

## Python implementation

The code below is a direct implementation of the equations in the Wikipedia article.

Note that longitude and latitude below are assumed to be in radians. You can convert from degrees to radians with SciPy’s `deg2rad`

function.

from scipy import sin, cos, tan, arctan, arctan2, arccos, pi def spherical_distance(lat1, long1, lat2, long2): phi1 = 0.5*pi - lat1 phi2 = 0.5*pi - lat2 r = 0.5*(6378137 + 6356752) # mean radius in meters t = sin(phi1)*sin(phi2)*cos(long1-long2) + cos(phi1)*cos(phi2) return r * arccos(t) def ellipsoidal_distance(lat1, long1, lat2, long2): a = 6378137.0 # equatorial radius in meters f = 1/298.257223563 # ellipsoid flattening b = (1 - f)*a tolerance = 1e-11 # to stop iteration phi1, phi2 = lat1, lat2 U1 = arctan((1-f)*tan(phi1)) U2 = arctan((1-f)*tan(phi2)) L1, L2 = long1, long2 L = L2 - L1 lambda_old = L + 0 while True: t = (cos(U2)*sin(lambda_old))**2 t += (cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda_old))**2 sin_sigma = t**0.5 cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(lambda_old) sigma = arctan2(sin_sigma, cos_sigma) sin_alpha = cos(U1)*cos(U2)*sin(lambda_old) / sin_sigma cos_sq_alpha = 1 - sin_alpha**2 cos_2sigma_m = cos_sigma - 2*sin(U1)*sin(U2)/cos_sq_alpha C = f*cos_sq_alpha*(4 + f*(4-3*cos_sq_alpha))/16 t = sigma + C*sin_sigma*(cos_2sigma_m + C*cos_sigma*(-1 + 2*cos_2sigma_m**2)) lambda_new = L + (1 - C)*f*sin_alpha*t if abs(lambda_new - lambda_old) <= tolerance: break else: lambda_old = lambda_new u2 = cos_sq_alpha*((a**2 - b**2)/b**2) A = 1 + (u2/16384)*(4096 + u2*(-768+u2*(320 - 175*u2))) B = (u2/1024)*(256 + u2*(-128 + u2*(74 - 47*u2))) t = cos_2sigma_m + 0.25*B*(cos_sigma*(-1 + 2*cos_2sigma_m**2)) t -= (B/6)*cos_2sigma_m*(-3 + 4*sin_sigma**2)*(-3 + 4*cos_2sigma_m**2) delta_sigma = B * sin_sigma * t s = b*A*(sigma - delta_sigma) return s

Reminds me of avionics work I did using the WMM (World Magnetic Model) to estimate the local magnetic heading based on GPS. We needed +/- 0.1 degree resolution with +/- .25 degree accuracy. The model itself was far too complex to run in real-time on our embedded host processor (https://www.ngdc.noaa.gov/geomag/WMM/back.shtml), but stacked LUTs and a fancy interpolator did the job just fine.

The testing and certification took many times longer than the algorithm development. The FAA inspector wanted to see simulated performance flying near the poles and the SAA (South Atlantic Anomaly). Many, many times. From every approach angle. Then in CW and CCW orbits. Ugh.

We also had tons of flight data from our electronic compass product, so we also replayed those. Numerous differences were seen, *all* of which were traced to errors in the electronic compass. (Mainly hard and soft iron errors.)

Wouldn’t it be fairer to the spherical estimate to use the mean radius rather than the equatorial? This should bring the underestimates and the overestimates closer to one another.

I wonder what value of radius of the spherical approximation, between shorter and longer radius of ellipsoid, would minimize errors…

Dr. Drang: I revised the post to use mean radius.

Jakub: I don’t know exactly, but I bet the optimal radius is not far from the mean radius. Maybe closer to the equatorial radius because that’s the radius in two orthogonal directions while the polar radius is the radius in one direction. Maybe the geometric mean, cube root of equatorial^2 * polar.