Fixing a Haversine Function That Returns the Wrong Distance

Problem statement

You wrote or copied a haversine function, and the distance between New York JFK and London Heathrow comes out as 7,681 km, or 8,186 km, or 0.87, instead of about 5,540 km. The formula is right; the inputs are not. Four mistakes account for almost every wrong haversine result: passing degrees to sin and cos, which expect radians; swapping latitude and longitude; using scikit-learn's haversine_distances, which returns distances on a sphere of radius 1; and mixing up the Earth's radius in kilometres and miles.

This page reproduces each mistake on real airport coordinates from OurAirports, shows the number it gives against the correct one, and adds the checks that catch them before the numbers reach a report.

Quick answer

Convert degrees to radians before any trigonometry, keep the argument order (lat1, lon1, lat2, lon2), and multiply angular distances by the radius in the unit you want:

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))

def haversine_degrees_bug(lat1, lon1, lat2, lon2, r=6371.0088):
    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))

j, l = AP.loc["JFK"], AP.loc["LHR"]
print(f"correct:            {haversine_km(j.latitude_deg, j.longitude_deg, l.latitude_deg, l.longitude_deg):,.1f} km")
print(f"degrees, no radians: {haversine_degrees_bug(j.latitude_deg, j.longitude_deg, l.latitude_deg, l.longitude_deg):,.1f} km")
print(f"exact (ellipsoid):  {geod.inv(j.longitude_deg, j.latitude_deg, l.longitude_deg, l.latitude_deg)[2] / 1000:,.1f} km")

correct:            5,539.8 km
degrees, no radians: 7,681.4 km
exact (ellipsoid):  5,554.7 km

With radians, JFK to Heathrow is 5,539.8 km, 0.3% short of the exact ellipsoidal 5,554.7 km. Without the conversion the same formula returns 7,681.4 km: a believable number, and 39% too long.

Degrees for radians, swapped order, unit-sphere results

1. Degrees instead of radians

Without the conversion, sin(51.47) treats 51.47 as radians, about eight full turns, and the result is essentially random. Measure how wrong it is across many pairs:

rng = np.random.default_rng(7)
i, k = rng.integers(0, len(AP), (2, 5000))
a, b = AP.iloc[i], AP.iloc[k]
good = haversine_km(a.latitude_deg.values, a.longitude_deg.values, b.latitude_deg.values, b.longitude_deg.values)
bad = haversine_degrees_bug(a.latitude_deg.values, a.longitude_deg.values, b.latitude_deg.values, b.longitude_deg.values)
ok = good > 1
print(f"pairs: {ok.sum():,}; correlation of buggy with correct distances: {np.corrcoef(good[ok], bad[ok])[0, 1]:.2f}")
print(f"buggy result within 10% of the truth: {(np.abs(bad[ok] / good[ok] - 1) < 0.1).mean():.1%}")

pairs: 4,994; correlation of buggy with correct distances: 0.02
buggy result within 10% of the truth: 9.2%

Across 4,994 random pairs of airports the buggy distances have a correlation of 0.02 with the correct ones: no relation at all. 9.2% happen to land within 10% of the truth, which is what makes a spot check unreliable.

2. Latitude and longitude swapped

Most data stores longitude first (x, y), most haversine functions take latitude first. Swapping gives distances that look plausible but are wrong:

sw = haversine_km(j.longitude_deg, j.latitude_deg, l.longitude_deg, l.latitude_deg)
print(f"JFK-LHR with lon and lat swapped: {sw:,.1f} km (correct {haversine_km(j.latitude_deg, j.longitude_deg, l.latitude_deg, l.longitude_deg):,.1f} km)")
bad_sw = haversine_km(a.longitude_deg.values, a.latitude_deg.values, b.longitude_deg.values, b.latitude_deg.values)
print(f"random pairs, swapped: median error {np.median(np.abs(bad_sw[ok] / good[ok] - 1)):.0%}; within 10%: {(np.abs(bad_sw[ok] / good[ok] - 1) < 0.1).mean():.0%}")
print(f"swapped latitudes outside -90..90: {(np.abs(AP.longitude_deg) > 90).mean():.0%} of airports")

JFK-LHR with lon and lat swapped: 8,185.8 km (correct 5,539.8 km)
random pairs, swapped: median error 17%; within 10%: 35%
swapped latitudes outside -90..90: 28% of airports

Swapping latitude and longitude gives 8,185.8 km for JFK to Heathrow. Over random pairs the swapped result is a median 17% off, and 35% of pairs are within 10%, close enough to pass a glance. One clue is in the data: 28% of airports have longitudes beyond ±90, which cannot be latitudes.

3. scikit-learn returns radians on a unit sphere

sklearn.metrics.pairwise.haversine_distances expects [latitude, longitude] in radians and returns the central angle. Multiply by the radius:

from sklearn.metrics.pairwise import haversine_distances
X = np.radians([[j.latitude_deg, j.longitude_deg], [l.latitude_deg, l.longitude_deg]])
d = haversine_distances(X)[0, 1]
print(f"sklearn output: {d:.6f} (an angle in radians); x 6,371.0088 = {d * 6371.0088:,.1f} km; x 3,958.76 = {d * 3958.7613:,.1f} miles")
Xdeg = np.array([[j.latitude_deg, j.longitude_deg], [l.latitude_deg, l.longitude_deg]])
print(f"same call with degrees: {haversine_distances(Xdeg)[0, 1] * 6371.0088:,.1f} km")

sklearn output: 0.869538 (an angle in radians); x 6,371.0088 = 5,539.8 km; x 3,958.76 = 3,442.3 miles
same call with degrees: 7,681.4 km

scikit-learn returns 0.869538, the central angle in radians. Times 6,371.0088 it is 5,539.8 km, the same as the correct function; times 3,958.76 it is 3,442.3 miles. Called with degrees, it returns the same wrong 7,681.4 km as the buggy function.

A test that catches all four mistakes

The data used on this page

Coordinates of the large airports with an IATA code in OurAirports, a public-domain database, and exact distances from pyproj.Geod on WGS84 as the reference.

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")

1,172 large airports

The data has 1,172 large airports with an IATA code.

Example 1: a test that catches all four

A function that passes these checks is very unlikely to have any of the bugs above:

ref = geod.inv(j.longitude_deg, j.latitude_deg, l.longitude_deg, l.latitude_deg)[2] / 1000
def check(f):
    checks = {
        "zero for the same point": abs(f(51.47, -0.45, 51.47, -0.45)) < 1e-9,
        "one degree of latitude ~111 km": abs(f(0, 0, 1, 0) - 111.2) < 0.5,
        "one degree of longitude at 60°N ~55.6 km": abs(f(60, 0, 60, 1) - 55.6) < 0.5,
        "JFK-LHR within 1% of the ellipsoid": abs(f(j.latitude_deg, j.longitude_deg, l.latitude_deg, l.longitude_deg) / ref - 1) < 0.01,
        "symmetric": abs(f(10, 20, 30, 40) - f(30, 40, 10, 20)) < 1e-9,
    }
    return {k: bool(v) for k, v in checks.items()}
print("correct function:", check(haversine_km))
print("degrees bug:", check(haversine_degrees_bug))
print("swapped arguments:", check(lambda a, b, c, d: haversine_km(b, a, d, c)))

correct function: {'zero for the same point': True, 'one degree of latitude ~111 km': True, 'one degree of longitude at 60°N ~55.6 km': True, 'JFK-LHR within 1% of the ellipsoid': True, 'symmetric': True}
degrees bug: {'zero for the same point': True, 'one degree of latitude ~111 km': False, 'one degree of longitude at 60°N ~55.6 km': False, 'JFK-LHR within 1% of the ellipsoid': False, 'symmetric': True}
swapped arguments: {'zero for the same point': True, 'one degree of latitude ~111 km': True, 'one degree of longitude at 60°N ~55.6 km': False, 'JFK-LHR within 1% of the ellipsoid': False, 'symmetric': True}

The correct function passes all five checks. The degrees bug passes only the two that any function of this shape passes, zero distance and symmetry, and fails every known distance. The swapped call passes the one-degree-of-latitude test by coincidence, because one degree of latitude at the equator and one degree of longitude at the equator are equal, but fails at 60°N and on the real route: test more than one known distance.

Trigonometry works in radians

NumPy's and Python's trigonometric functions work in radians, a full turn being 2π, about 6.28. A latitude of 51.47 passed straight to np.sin is read as 51.47 radians, a little over eight full turns, so the sine has nothing to do with the latitude, and the result is a number of the right size, thousands of kilometres, with no relation to the real distance. That is what makes the bug dangerous: it produces plausible-looking output. The same happens with a swap of latitude and longitude, except that the result is often close to right for short distances and badly wrong for long ones, because the formula weighs the two coordinates differently.

The unit-sphere convention of scikit-learn is deliberate, so its BallTree and pairwise functions work for any sphere; it leaves the multiplication by the radius to you. The radius itself is the last trap: 6,371 for kilometres, 3,959 for miles, 3,440 for nautical miles, and mixing them gives errors of a factor of 1.6 or 1.85 that a sanity check catches at once.

math against numpy and known-distance tests

• math vs numpy: math.radians works on single numbers, np.radians on arrays; both are fine.
• Pandas columns: pass .values or Series of the same index; misaligned indexes silently pair the wrong rows.
• Antipodal rounding: for nearly opposite points, sqrt(a) can slightly exceed 1 through rounding; clip it to 1 before arcsin.
• Known-distance tests: one degree of latitude is about 111.2 km everywhere; one degree of longitude is 111.2 × cos(latitude).

Internal links

• How to calculate the distance between two coordinates in Python — the correct version in full
• How to find all points within a radius in Python — scikit-learn's unit sphere again
• Latitude and longitude explained — which number is which

FAQ

Why does my haversine function give the wrong distance?

Usually because the coordinates were not converted to radians, latitude and longitude are swapped, or a unit-sphere result was not multiplied by the Earth's radius.

Does np.sin take degrees?

No, radians; convert with np.radians first.

Why does sklearn's haversine_distances return a number below 3.15?

It returns the central angle in radians on a unit sphere; multiply by 6,371.0088 for kilometres.

How can I test a distance function?

Check that one degree of latitude is about 111.2 km, a known route matches a reference to within 1%, the distance from a point to itself is zero and the function is symmetric.