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.

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.

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.

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.