How to Calculate the Distance Between Two Coordinates in Python (Haversine)
On this page
Problem statement
Two points given as latitude and longitude, and you need the distance between them in kilometres: from a customer to a store, between two airports, along each step of a GPS track. Subtracting the coordinates and applying Pythagoras gives nonsense, because degrees of longitude shrink towards the poles and the Earth is round. The haversine formula gives the great-circle distance on a sphere in a few lines of NumPy, and it is accurate to about half a percent either way, which is enough for most purposes once you know where that error comes from.
This page implements the haversine formula, applies it to real airport coordinates from OurAirports one pair at a time and for thousands of pairs at once, measures its error against the exact distance on the WGS84 ellipsoid, and shows what the choice of Earth radius does.
Quick answer
Convert degrees to radians, apply the haversine formula, multiply by the Earth's mean radius:
import numpy as np
def haversine_km(lat1, lon1, lat2, lon2, r=6371.0088):
lat1, lon1, lat2, lon2 = map(np.radians, (lat1, lon1, lat2, lon2))
a = np.sin((lat2 - lat1) / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin((lon2 - lon1) / 2) ** 2
return 2 * r * np.arcsin(np.sqrt(a))
jfk, lhr = AP.loc["JFK"], AP.loc["LHR"]
print(f"JFK to LHR: {haversine_km(jfk.latitude_deg, jfk.longitude_deg, lhr.latitude_deg, lhr.longitude_deg):,.1f} km")JFK to LHR: 5,539.8 km
New York JFK to London Heathrow is 5,539.8 km along the great circle.
Many distances at once, accuracy, which radius
1. Many distances at once
The function works on NumPy arrays, so distances from one point to every airport, or between two columns, need no loop:
d = haversine_km(lhr.latitude_deg, lhr.longitude_deg, AP.latitude_deg.values, AP.longitude_deg.values)
near = pd.Series(d, index=AP.index).drop("LHR").nsmallest(5)
print("large airports nearest to Heathrow:", ", ".join(f"{k} {v:.0f} km" for k, v in near.items()))
print(f"distances computed: {len(d):,}; farthest from LHR: {pd.Series(d, index=AP.index).idxmax()} at {d.max():,.0f} km")large airports nearest to Heathrow: LGW 41 km, LTN 45 km, STN 66 km, BHX 140 km, BRS 157 km
distances computed: 1,172; farthest from LHR: CHC at 18,984 km
One call computes the distance from Heathrow to all 1,172 airports. The nearest large airports are Gatwick, 41 km away, Luton, 45 km, and Stansted, 66 km; the farthest is Christchurch, New Zealand, 18,984 km away.
2. How accurate is it?
The haversine formula assumes a sphere; the Earth is an ellipsoid, flattened by about 1 part in 298. Compare with the exact geodesic distance on WGS84 for 20,000 random pairs of airports:
rng = np.random.default_rng(1)
i, j = rng.integers(0, len(AP), (2, 20000))
a, b = AP.iloc[i], AP.iloc[j]
hav = haversine_km(a.latitude_deg.values, a.longitude_deg.values, b.latitude_deg.values, b.longitude_deg.values)
geo = geod.inv(a.longitude_deg.values, a.latitude_deg.values, b.longitude_deg.values, b.latitude_deg.values)[2] / 1000
ok = geo > 1
err = (hav[ok] - geo[ok]) / geo[ok] * 100
print(f"pairs: {ok.sum():,}; error against the ellipsoid: median {np.median(err):+.3f}%, range {err.min():+.3f}% to {err.max():+.3f}%")
print(f"largest absolute difference: {np.abs(hav[ok] - geo[ok]).max():.1f} km")pairs: 19,981; error against the ellipsoid: median -0.095%, range -0.379% to +0.555%
largest absolute difference: 37.6 km
Against the exact distance on the WGS84 ellipsoid, the haversine distance is 0.095% short for the median pair, and the errors span -0.379% to +0.555%. The largest difference in kilometres is 37.6 km, on a long route. The sign depends on where and in which direction a route runs: near the equator, the mean sphere is too small for east-west arcs and too large for north-south ones.
3. Which radius?
The Earth has no single radius: 6,378.137 km at the equator, 6,356.752 km at the poles. The mean radius, 6,371.0088 km, minimises the average error. Compare the choices on one long route:
sin = AP.loc["SIN"]
exact = geod.inv(jfk.longitude_deg, jfk.latitude_deg, sin.longitude_deg, sin.latitude_deg)[2] / 1000
for name, r in (("mean (IUGG)", 6371.0088), ("equatorial", 6378.137), ("polar", 6356.752), ("rounded 6,371", 6371.0)):
h = haversine_km(jfk.latitude_deg, jfk.longitude_deg, sin.latitude_deg, sin.longitude_deg, r)
print(f"{name:>14}: JFK to SIN {h:,.1f} km ({h - exact:+.1f} km against the ellipsoid, {exact:,.1f} km)") mean (IUGG): JFK to SIN 15,340.6 km (-8.1 km against the ellipsoid, 15,348.7 km)
equatorial: JFK to SIN 15,357.8 km (+9.1 km against the ellipsoid, 15,348.7 km)
polar: JFK to SIN 15,306.3 km (-42.4 km against the ellipsoid, 15,348.7 km)
rounded 6,371: JFK to SIN 15,340.6 km (-8.1 km against the ellipsoid, 15,348.7 km)
On New York to Singapore, 15,348.7 km on the ellipsoid, the mean radius is 8.1 km short, the equatorial radius 9.1 km long and the polar radius 42.4 km short. Rounding the mean radius to 6,371 km changes nothing at this precision.
A distance matrix of airports and the farthest pair
The data used on this page
Airport coordinates from OurAirports, a public-domain database, keeping the large airports that have an IATA code. The exact ellipsoidal distances come from pyproj.Geod on WGS84.
from pathlib import Path
import numpy as np
import pandas as pd
import requests
from pyproj import Geod
Path("data").mkdir(exist_ok=True)
p = Path("data/airports.csv")
if not p.exists():
r = requests.get("https://davidmegginson.github.io/ourairports-data/airports.csv", timeout=300)
r.raise_for_status()
p.write_bytes(r.content)
AP = pd.read_csv(p, keep_default_na=False, na_values=[""])
AP = AP[(AP.type == "large_airport") & AP.iata_code.notna()].set_index("iata_code")
geod = Geod(ellps="WGS84")
print(f"{len(AP):,} large airports with an IATA code; JFK at {AP.loc['JFK', 'latitude_deg']:.4f}, {AP.loc['JFK', 'longitude_deg']:.4f}")1,172 large airports with an IATA code; JFK at 40.6394, -73.7793
OurAirports lists 1,172 large airports with an IATA code; JFK sits at 40.6394°N, 73.7793°W.
Example 1: a distance matrix and the farthest pair
For a few hundred points, a full matrix by broadcasting is fast. Find the two large airports farthest apart:
lat, lon = AP.latitude_deg.values, AP.longitude_deg.values
M = haversine_km(lat[:, None], lon[:, None], lat[None, :], lon[None, :])
i, j = np.unravel_index(M.argmax(), M.shape)
print(f"matrix {M.shape[0]} x {M.shape[1]}; farthest pair: {AP.index[i]} ({AP.municipality.iloc[i]}) and {AP.index[j]} ({AP.municipality.iloc[j]}), {M[i, j]:,.0f} km")
print(f"half the Earth's circumference: {np.pi * 6371.0088:,.0f} km")matrix 1172 x 1172; farthest pair: ESM (Tachina) and PDG (Padang (Katapiang)), 19,991 km
half the Earth's circumference: 20,015 km
The 1,172 × 1,172 matrix, 1.4 million distances, takes one broadcast. The farthest pair is Esmeraldas in Ecuador and Padang in Indonesia, 19,991 km apart, within 24 km of half the Earth's circumference: the two airports are nearly antipodal.
What the haversine formula computes
The haversine formula computes the central angle between two points on a sphere from their latitudes and longitudes; multiplied by the radius it gives the length of the great-circle arc, the shortest path along the surface. It is preferred to the spherical law of cosines because it stays accurate for points very close together, where the cosine form loses precision to rounding. The name comes from the haversine function, hav(θ) = sin²(θ/2), which made the formula convenient with printed tables.
Its limit is the sphere. The Earth's equatorial radius is 21 km longer than its polar radius, so no single radius fits every direction: a north-south arc and an east-west arc of the same angle differ in length. With the mean radius, the haversine distances between 19,981 random pairs of airports had a median error of -0.10% against the exact ellipsoidal value, and the extremes were -0.38% and +0.56%. For navigation, surveying or legal boundaries, use a geodesic on the ellipsoid, as pyproj.Geod or geopy.distance.geodesic compute; for analytics, nearest-neighbour searches and anything feeding a model, haversine is fast and good enough.
Radians, argument order and road distance
- Degrees, not radians: forgetting the conversion is the most common bug; see the radians fix.
- Order of arguments: latitude first or longitude first; mixing them gives plausible but wrong distances.
- Altitude: haversine ignores height; for aircraft at 11 km the difference in path length is under 0.2%.
- Road distance: great-circle distance is a lower bound for travel; use a routing engine for driving distance.
Internal links
- How to measure distance accurately in Python — the exact ellipsoidal distance
- How to find points within a radius in Python — haversine for searches
- Great circle and rhumb line explained — which path the distance measures
FAQ
How do I calculate the distance between two latitude/longitude points in Python?
Convert the coordinates to radians and use the haversine formula: a = sin²(Δφ/2) + cos φ1 cos φ2 sin²(Δλ/2), distance = 2R·arcsin(√a), with R = 6,371 km.
How accurate is the haversine formula?
Within about 0.6% of the exact distance on the WGS84 ellipsoid; for 19,981 random pairs of airports the median error was -0.10% and the worst +0.56%, 37.6 km on a long route.
What Earth radius should I use?
The mean radius, 6,371.0088 km (often rounded to 6,371 km); use 3,958.8 miles for miles.
Is there a Python library for haversine?
Yes: the haversine package, sklearn.metrics.pairwise.haversine_distances (in radians, on the unit sphere) and, for exact distances, pyproj.Geod or geopy.distance.