warming up your workspace

The haversine formula, or why you cannot just use Pythagoras on a map

You have two GPS coordinates and you want the distance between them. The obvious move is to treat latitude and longitude like x and y, subtract, and run Pythagoras. For two points across a city that is nearly fine. For two points across a continent it can be hundreds of kilometers wrong, and worse, it is wrong in a way that changes with where you are on the planet. The reason is simple: latitude and longitude are angles on a sphere, not meters on a plane. The formula that fixes it is called haversine, and it is short.

The one idea

A degree of latitude is always about 111 kilometers. A degree of longitude is not. At the equator it is also 111 kilometers, but near the poles the meridians crowd together and a degree of longitude shrinks toward zero. Pythagoras assumes both are constant, so it overcounts east-west distance everywhere except the equator. Haversine works on the actual sphere: it measures the angle between two points as seen from the center of the Earth, then multiplies by the radius to get arc length along the surface. That arc is the great-circle distance, the shortest path a plane can fly.

Build it

The formula takes the two latitudes and the differences in latitude and longitude, all in radians, and returns the surface distance.

from math import radians, sin, cos, asin, sqrt

def haversine(lat1, lon1, lat2, lon2, R=6371.0088):   # mean Earth radius, km
    phi1, phi2 = radians(lat1), radians(lat2)
    dphi = radians(lat2 - lat1)
    dlam = radians(lon2 - lon1)
    a = sin(dphi/2)**2 + cos(phi1)*cos(phi2)*sin(dlam/2)**2
    return 2 * R * asin(sqrt(a))

Three details that matter:

  • The cos(phi1) * cos(phi2) term is the whole correction. It shrinks the longitude contribution as you move away from the equator, exactly the effect Pythagoras misses. Take it out and you are back to a flat map.
  • Everything happens in radians. Latitude and longitude arrive in degrees, so the first thing the function does is convert. Forgetting this is the single most common bug in distance code.
  • The name comes from the "haversine" function, sin squared of half the angle. It is written this way instead of a plain cosine because for two very close points the plain-cosine version loses precision, the numbers round to the same value, while this form stays accurate down to meters.

Proof: it matches reality, and the flat map does not

jfk = (40.6413, -73.7781)     # New York
lhr = (51.4700, -0.4543)      # London
d = haversine(*jfk, *lhr)
print(f"JFK to LHR: {d:.0f} km")
print(f"quarter turn on equator: {haversine(0,0, 0,90):.0f} km")

Running it prints:

JFK to LHR: 5540 km
quarter turn on equator: 10008 km

The New York to London distance matches what an airline publishes. The equator check is a proof by geometry: a quarter of the way around the planet should be a quarter of the circumference, and 10,008 kilometers is exactly that. Now compare the flat-map shortcut on the same transatlantic route: it returns 5,785 kilometers, off by 245 kilometers, roughly the width of England, because it counted the longitude gap as if New York and London sat on the equator. For a delivery-radius check or a "nearest store" query, that error is the difference between a right answer and a wrong one.

Where this shows up

Every "find restaurants near me", every ride-hailing fare estimate, every geofence, and every shipping-zone calculation is doing this under the hood, or should be. Haversine assumes a perfect sphere, which is accurate to about 0.5 percent; when you need survey-grade precision you move up to Vincenty's formula on an ellipsoid, but the idea is the same, measure the angle, then walk the curve.

If you want to build the rest of what sits on top of coordinates, projections, tiles, spatial indexes, and the map math that turns a GPS pin into a place, that is what the geospatial track on IWTLP constructs from the ground up.