 johndcook.com Go back Open original

# Ellipsoid distance on Earth

By John 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%).

## 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 a = 6378137.0 # equatorial radius in meters
f = 1/298.257223563 # ellipsoid flattening
b = (1 - f)*a
tolerance = 1e-11 def spherical_distance(lat1, long1, lat2, long2): phi1 = 0.5*pi - lat1 phi2 = 0.5*pi - lat2 t = sin(phi1)*sin(phi2)*cos(long1-long2) + cos(phi1)*cos(phi2) return a * arccos(t) def ellipsoidal_distance(lat1, long1, lat2, long2): 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
```