Haversine Formula

The goal for this formula is to calculate the shortest great circle distance between two points on the globe designated by latitude and longitudes. The added benefit of the Haversine equation is that it calculates the central angle as well where $s = r\theta$.

The Haversine formula assumes a constant radius Earth and does not account for elevation changes between locations.

$$ d = 2r \arcsin\left( \sqrt{\sin^2 \left(\frac{\phi_2-\phi_1}{2}\right) + \cos(\phi_1)\cos(\phi_2)\sin^2 \left( \frac{\lambda_2-\lambda_1}{2} \right) } \right) $$

where:

  • $\phi$: latitude
  • $\lambda$: longitude
  • $r$: Radius of the sphere

In [1]:
import numpy as np
import math

# Mean radius of the earth
EARTH_RADIUS = 6371.009

def haversine(lat1, lon1, lat2, lon2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    
    Return (central angle, distance between points in km)
    """
    # convert decimal degrees to radians
    lat1, lon1, lat2, lon2 = [math.radians(x) for x in [lat1, lon1, lat2, lon2]]

    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    
    a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    central_angle = 2 * math.asin(math.sqrt(a))

    # s = r * theta
    km = EARTH_RADIUS * central_angle
    return (central_angle, km)

Test around the equator

Test the Haversine formula against a half distance around the equator. Assumes the equitorial radius is equal to the mean radius.


In [2]:
print(haversine(0, -90, 0, 90))
print(2 * np.pi * EARTH_RADIUS / 2)


(3.141592653589793, 20015.115070354455)
20015.115070354455

Test Against the Great Circle Mapper

  • From: DXB 25°15'10"N 55°21'52"E
  • To: SFO 37°37'08"N 122°22'32"W
  • Distance: 13041 km
The Great Circle Mapper uses the same WGS 84 reference ellipsoid used by the Global Positioning System (GPS). The WGS 84 ellipsoid has a polar radius of 6356.752 km and an equatorial radius of 6378.137 km.

ref: http://www.gcmap.com/faq/gccalc#ellipsoid


In [3]:
lat1 = 25+15/60.+10/3600.
lon1 = 55+21/60.+52/3600.
lat2 = 37+37/60.+8/3600.
lon2 = -122+22/60.+32/3600.
haversine(lat1, lon1, lat2, lon2)


Out[3]:
(2.0431640443445644, 13017.01651499562)

In [3]: