Why a Distance Calculation From GPS Coordinates Was Wrong by 20%

#ai #debugging #llm #machinelearning

Your distances are consistently too long, the error is roughly proportional, and it gets worse the further from the equator you test. That combination has a small number of causes, and the size of the error tells you which one.

The symptom, and how to triage it in one minute

Distances computed from latitude and longitude disagree with a mapping service, a GPS device or a known ground truth. Triage on the shape of the disagreement:

Off by 10–60%, worse at higher latitude, worse on east–west legs — Pythagoras on degrees. Section three.
Off by a clean factor of roughly 1/cos(latitude), on every leg equally — you are measuring in Web Mercator metres. Section four.
Off by well under 1% — sphere versus ellipsoid, and it may not be a bug at all. Section five.
Off by thousands of kilometres, or your point is in the Gulf of Guinea — latitude and longitude are swapped, or a missing value became 0. GeoJSON (RFC 7946) orders a position as longitude, latitude; most consumer APIs and most people say latitude first.

Reproducing the 20%

Two points in London, chosen so the leg runs diagonally:

A = 51.5074 N,  -0.1278 E
B = 51.5374 N,  -0.0978 E      (dLat = 0.0300 deg, dLon = 0.0300 deg)

WRONG: Pythagoras on degrees, then scale by 111,320 m/deg
  sqrt(0.0300^2 + 0.0300^2)  = 0.0424264 deg
  0.0424264 * 111320         = 4,723 m

RIGHT: haversine, R = 6,371,008.8 m
  a = sin^2(dLat/2) + cos(latA)cos(latB) sin^2(dLon/2)
    = 6.854e-8 + (0.62239 * 0.62198 * 6.854e-8)
    = 6.854e-8 + 2.653e-8  =  9.507e-8
  d = 2 R asin(sqrt(a))     = 3,929 m

error  4,723 / 3,929 = 1.202   ->  20.2% too long

Twenty per cent, from a coordinate pair you can paste into anything to check. The exact figure is not universal — it depends on latitude and on the bearing of the leg — and that dependence is the whole explanation.

Cause 1: degrees are not metres

A degree of latitude is about 111,320 m everywhere, because meridians are all great circles. A degree of longitude is 111,320 × cos(latitude) metres, because parallels shrink toward the poles. At 51.5074° that cosine is 0.62239, so a degree of longitude is 69,290 m — 62% of a degree of latitude.

Applying one scale factor to both axes therefore inflates the east–west component by 1 / cos(latitude) while leaving the north–south component correct. The total error depends on the bearing:

error factor = sqrt(dLat^2 + dLon^2) / sqrt(dLat^2 + (dLon * cos(lat))^2)

at latitude 51.5074 (cos = 0.62239):
  due north   (dLon = 0)          ->  1.000     0% error
  45 deg      (dLat = dLon)       ->  1.201    20% error   <- the worked case
  due east    (dLat = 0)          ->  1.607    61% error

Which explains the reports. Someone testing in Singapore at 1.3°N sees a 0.03% error and concludes the code is fine; the same code in Oslo at 59.9° is 100% wrong on an east–west leg. And an aggregate over many legs lands somewhere in between, which is why the total looks like a mysterious constant fudge factor rather than a bug.

There is one more twist worth knowing, because it produces this exact 20% from a different cause. Swap latitude and longitude and feed the pair to a correct haversine: the latitude terms become −0.1278 and −0.0978, both near the equator, where cos(lat) ≈ 1.0. The cosine term stops doing anything, and haversine returns 4,723 m — the same number as the broken planar formula. A correct function with transposed inputs and a broken function with correct inputs give identical answers here, so a matching result proves nothing. Check a north–south leg too, where the two diverge.

Cause 2: metres that are not metres

The other common version: the coordinates are projected, so the units really are metres, and the distance is still wrong. If the projection is Web Mercator, EPSG:3857 — the default of most web mapping stacks, and what PostGIS gives you if you reproject “to metres” without thinking — then distances are inflated by the local scale factor, which is 1 / cos(latitude).

true ground distance                              3,929 m
Web Mercator scale factor at 51.5074  = 1/0.62239 = 1.607
distance measured in EPSG:3857 metres             = 6,313 m   61% too long

Mercator is conformal, so that scale factor is the same in every direction at a point — which makes the error a clean multiplier rather than a bearing-dependent one, and that is the diagnostic that separates this from cause 1. It is also why Greenland looks the size of Africa. Web Mercator is a display projection; it is not a measurement projection, and there is no configuration that makes it one.

The correct projected choice is a local one: the appropriate UTM zone, a national grid such as EPSG:27700 for Great Britain, or a state-plane zone. Those are designed so the scale factor stays within a small fraction of a per cent across the zone. Cross a zone boundary and you are back to needing geodesics.

Cause 3: sphere versus ellipsoid

Haversine is not the last word; it assumes a sphere. The WGS84 ellipsoid that GPS reports positions on has a semi-major axis of 6,378,137.0 m and a flattening of 1/298.257223563, giving a semi-minor axis of about 6,356,752.3 m — the poles are roughly 21.4 km closer to the centre than the equator, about 0.34% of the mean radius.

So haversine differs from a true WGS84 geodesic by up to roughly 0.5%, depending on latitude and bearing. On the 3,929 m leg above that is on the order of 20 m. Whether that matters is a requirements question, not a correctness one: it is invisible in a delivery-zone area calculation and unacceptable in a survey or an aviation leg. Where it does matter, use an inverse geodesic on the ellipsoid — Karney’s algorithm, implemented in GeographicLib, and exposed as ST_Distance on the PostGIS geography type and as Geod.inv in pyproj — which is accurate to well under a millimetre and has none of Vincenty’s convergence failures on near-antipodal pairs.

The fix

Decide the tolerance first. Under 0.5% error acceptable and legs under a few hundred kilometres: haversine is enough. Tighter than that, or antipodal-ish legs: use a geodesic.
Never apply a planar formula to unprojected degrees. If the data is in EPSG:4326, either use a spherical or ellipsoidal function, or project to a local metric CRS first — and confirm the target CRS is not EPSG:3857.
In PostGIS, store the column as geography rather than geometry when the answers should be in metres on the globe. ST_Distance on geography is a spheroidal distance in metres; on geometry with SRID 4326 it is a number of degrees that looks like a distance and is not one.
Pin coordinate order at every boundary. Write down whether your internal tuple is (lat, lon) or (lon, lat), name the fields rather than using positional tuples where you can, and assert that latitude is within ±90 on ingest — the single cheapest check, since it catches every swap outside the tropics.
Regression-test with two legs at a mid-to-high latitude: one north–south and one east–west, with known answers. One diagonal leg alone will not distinguish the failures on this page, and neither will any test at the equator. The same pair of legs is worth reusing for geofence radius checks and for side-of-road offsets, which fail for exactly the same reason.