Detecting Geocoding Outliers with Haversine Distance

TL;DR: Catch plausible-but-wrong geocodes by measuring the haversine (great-circle) distance between the returned point and an expected reference — a postal centroid or the prior stored value — and flagging anything beyond an area-type threshold. This is the geometric semantic check behind validating geocoding accuracy and confidence scoring.

Why Distance Catches What Range Checks Miss

A coordinate can be individually valid on both axes and still be wrong. A result that lands in the correct country but 60 km from the postal code it claims is a classic provider error: it passes every range and bounding-box check yet is useless for delivery. The only way to catch it is to compare the result against an independent expectation and measure how far apart they are. Haversine distance gives that measurement over a sphere, needs no external libraries, and is fast enough to run over an entire batch.

The reference point is whatever you trust more than the new geocode: a centroid derived from the record’s postal code, the previously stored coordinate for that address, or a rooftop result from a higher-tier provider. If the new result strays too far from that anchor, flag it.

What a Distance Check Catches That a Range Check Cannot

Bounds checks confirm a coordinate is on the planet; they say nothing about whether it is in the right place. Three common failures produce perfectly legal latitude and longitude values, and only a distance comparison against an expected reference point exposes them.

The middle case is worth a dedicated assertion of its own, because it is both the most common and the easiest to miss: exactly zero is a legal coordinate and an overwhelmingly likely default. Treat (0, 0) as null unless your data genuinely includes the Gulf of Guinea.

Formula Terms and Edge Cases

The haversine formula computes the great-circle distance d between two latitude/longitude points on a sphere of radius R:

Two edge cases dominate in practice. The antimeridian: points either side of the ±180° meridian have longitudes that look far apart numerically, but the atan2-based c computes the true short-way angular distance correctly, so no manual unwrapping is needed. Missing reference: many records have no postal centroid or prior value; treat a missing reference as “cannot evaluate” and skip, never as an automatic failure.

Pure-Python Haversine

No dependencies — just the standard library math module. The function returns kilometres and clamps a to guard against floating-point drift pushing it slightly above 1.0 for antipodal points.

from __future__ import annotations

import math
from typing import Optional, Tuple

_EARTH_RADIUS_KM = 6371.0088

Coord = Tuple[float, float]  # (lat, lon) in decimal degrees


def haversine_km(a: Coord, b: Coord) -> float:
    """Great-circle distance between two (lat, lon) points, in kilometres.

    Args:
        a: First point as (latitude, longitude) in decimal degrees.
        b: Second point as (latitude, longitude) in decimal degrees.

    Returns:
        Distance in kilometres.

    Raises:
        ValueError: If either latitude or longitude is out of range.
    """
    lat1, lon1 = a
    lat2, lon2 = b
    for lat in (lat1, lat2):
        if not -90.0 <= lat <= 90.0:
            raise ValueError(f"latitude out of range: {lat}")
    for lon in (lon1, lon2):
        if not -180.0 <= lon <= 180.0:
            raise ValueError(f"longitude out of range: {lon}")

    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    d_phi = math.radians(lat2 - lat1)
    d_lambda = math.radians(lon2 - lon1)

    h = (
        math.sin(d_phi / 2) ** 2
        + math.cos(phi1) * math.cos(phi2) * math.sin(d_lambda / 2) ** 2
    )
    h = min(1.0, max(0.0, h))  # guard against FP drift
    c = 2 * math.atan2(math.sqrt(h), math.sqrt(1 - h))
    return _EARTH_RADIUS_KM * c


def is_outlier(
    result: Coord,
    reference: Optional[Coord],
    threshold_km: float,
) -> Optional[bool]:
    """Flag a geocode as an outlier if it lies beyond threshold_km of reference.

    Args:
        result:       Geocoded (lat, lon) to test.
        reference:    Trusted (lat, lon), e.g. a postal centroid, or None.
        threshold_km: Maximum acceptable distance for this area type.

    Returns:
        True if the result is an outlier, False if within tolerance, or
        None if no reference is available to evaluate against.
    """
    if reference is None:
        return None
    return haversine_km(result, reference) > threshold_km

Usage against a postal centroid:

# Provider returned a point; postal centroid is the trusted anchor.
result = (48.8600, 2.3400)
postal_centroid = (48.8566, 2.3522)
flagged = is_outlier(result, postal_centroid, threshold_km=5.0)
print(haversine_km(result, postal_centroid))  # ~0.93 km -> within tolerance
print(flagged)                                # False

Choosing a Threshold by Area Type

A single global threshold either misses urban errors or floods rural records with false positives. Postal codes cover vastly different areas by density, so scale the threshold to the area type.

When the reference is a prior stored coordinate rather than a centroid, tighten the threshold sharply — a genuine address does not move, so even a 500 m jump between runs signals a provider change worth reviewing.

Thresholds Follow Density, Not Taste

A single distance threshold is wrong everywhere. Two kilometres from a postcode centroid is a serious anomaly in Manhattan and entirely normal in rural Montana, where a single postcode can span a hundred kilometres. Tie the threshold to the area type the record falls in.

Deriving the area type is easier than it sounds: a population-density lookup or even the count of addresses sharing the postcode in your own data is enough to bucket records into these four classes. A fixed global threshold set for the urban case buries the rural team in false positives, and set for the rural case misses genuine urban errors entirely.

numpy-Vectorized pandas Variant

For batch screening, vectorize haversine over two coordinate columns and a reference pair. The array form computes millions of distances per second.

import numpy as np
import pandas as pd

_EARTH_RADIUS_KM = 6371.0088


def haversine_series(
    lat: pd.Series,
    lon: pd.Series,
    ref_lat: pd.Series,
    ref_lon: pd.Series,
) -> pd.Series:
    """Vectorized great-circle distance (km) between result and reference.

    NaN in any reference cell yields NaN distance (unevaluable), so records
    without a reference are naturally excluded from outlier flagging.
    """
    lat1 = np.radians(lat.to_numpy(dtype="float64"))
    lat2 = np.radians(ref_lat.to_numpy(dtype="float64"))
    d_phi = lat2 - lat1
    d_lambda = np.radians(ref_lon.to_numpy(dtype="float64")
                          - lon.to_numpy(dtype="float64"))
    h = (np.sin(d_phi / 2.0) ** 2
         + np.cos(lat1) * np.cos(lat2) * np.sin(d_lambda / 2.0) ** 2)
    h = np.clip(h, 0.0, 1.0)
    c = 2.0 * np.arctan2(np.sqrt(h), np.sqrt(1.0 - h))
    return pd.Series(_EARTH_RADIUS_KM * c, index=lat.index)


def flag_outliers(
    df: pd.DataFrame,
    threshold_col: str = "threshold_km",
) -> pd.DataFrame:
    """Add 'ref_distance_km' and 'is_outlier' columns to a results DataFrame.

    Expects columns: lat, lon, ref_lat, ref_lon, and a per-row threshold_km.
    Rows with a missing reference get is_outlier = False (unevaluable).
    """
    out = df.copy()
    dist = haversine_series(out["lat"], out["lon"],
                            out["ref_lat"], out["ref_lon"])
    out["ref_distance_km"] = dist.round(3)
    out["is_outlier"] = (dist > out[threshold_col]).fillna(False)
    return out

Carrying a per-row threshold_km column lets a single pass apply urban, suburban, and rural tolerances at once — derive it from the record’s area type during parsing rather than hard-coding one global value.

Vectorising Without Losing Precision

Haversine on a million rows is a numpy expression, not a loop — but two details decide whether the vectorised version is also correct. Radians must be converted once rather than per term, and the arcsine form is numerically better behaved than the naive arccosine for the small distances that dominate address work.

The precision note is not academic here. Urban thresholds live in the hundreds of metres, and an implementation that loses several digits at short range will happily report a forty-metre discrepancy as zero — hiding precisely the small-scale disagreements the check was added to surface.

Edge Cases

Antimeridian-spanning pairs

A result at 179.9°E and a reference at 179.9°W are about 22 km apart, not most of the way around the globe. The atan2 form used above returns the correct short-way distance without any manual longitude unwrapping, so Pacific and Fiji-area records need no special handling.

No reference available

Records with no postal centroid and no prior value cannot be distance-checked. The pure-Python is_outlier returns None and the vectorized variant leaves is_outlier False — never invent a reference, because a wrong anchor produces worse false positives than skipping.

Integration Note

Distance-to-reference is the sharpest form of the semantic tier in validating geocoding accuracy and confidence scoring: a result can carry a high composite from calculating geocoding confidence scores in Python and still be an outlier that lands kilometres from its postal centroid. Run both checks and require a result to clear each: a strong composite score AND a small reference distance before you commit the geocode.

---

Related

• Validating geocoding accuracy and confidence scoring — the three-tier validator where distance screening is the semantic check.
• Calculating geocoding confidence scores in Python — the composite score that pairs with this geometric test.
• Building regression corpora for geocoding accuracy — capture flagged outliers as fixtures so the same errors are caught in CI.