Skip to content
scipython.com

P4.4.2The haversine formula

The Haversine formula gives the shortest (great-circle) distance, \(d\), between two points on a sphere of radius \(R\) from their longitudes \((\lambda_1, \lambda_2)\) and latitudes \((\phi_1, \phi_2)\):

\[d = 2R \arcsin\left( \sqrt{\mathrm{haversin}(\phi_2-\phi_1) + \cos\phi_1\cos\phi_2\mathrm{haversin}(\lambda_2-\lambda_1)} \right),\]

where the haversine function of an angle is defined by

\[\mathrm{haversin}(\alpha) = \sin^2\left(\frac{\alpha}{2}\right).\]

Write a program to calculate the shortest distance in km between two points on the surface of the Earth (considered as a a sphere of radius 6378.1 km) given as two command line arguments, each of which is a comma-separated pair of latitude, longitude values in degrees. For example, the distance between Paris (48.9 °N, 2.4 °E)and Rome (41.9 °N, 12.5 °E) is given by executing:

$ python greatcircle.py 48.9,2.4 41.9,12.5
1107 km
Show solution

The solution below is one approach. Note that we must convert the provided arguments from string values representing angles in degrees to floats in radians.

import sys
import math

# Calculate the shortest distance between two points on a sphere provided on
# the command line as (longitude, latitude pairs). For example, Paris and Rome:
# 48.9,2.4 41.9,12.5

haversin = lambda alpha: math.sin(alpha/2)**2

def gc_distance(loc1, loc2, R=6378.1):
    """
    Return the great-circle distance between two points on a sphere of radius
    R, provided as (latitude, longitude) pairs, loc=(phi, lamdbda) in radians.

    """

    (phi1, lambda1), (phi2, lambda2) = loc1, loc2
    d = 2 * R * math.asin(math.sqrt(haversin(phi2-phi1)
                + math.cos(phi1)*math.cos(phi2)*haversin(lambda2-lambda1)))
    return d

loc1, loc2 = sys.argv[1:]
phi1, lambda1 = [math.radians(float(x)) for x in loc1.split(',')]
phi2, lambda2 = [math.radians(float(x)) for x in loc2.split(',')]

d = gc_distance((phi1, lambda1), (phi2, lambda2))
print('{} km'.format(int(d)))