Detecting Geocoding Outliers with Spatial Clustering

This page covers one validation task: finding the geocodes in a store or customer table that are almost certainly wrong — transposed coordinates, ZIP-centroid fallbacks, or points in the wrong hemisphere — using spatial density rather than eyeballing a map.

Coordinate quality is the foundation the rest of the location intelligence architecture and data foundations stands on — a wrong point poisons every join built on it. Geocoders fail quietly. A batch of a few thousand addresses will contain a handful where the provider returned the centroid of a ZIP code instead of a rooftop, where a downstream export swapped latitude and longitude, or where a null island coordinate (0, 0) slipped through. None of these throw an error; they produce a valid-looking point in the wrong place, and that point silently corrupts every trade area and demographic join built on top of it. Spatial clustering catches them because a genuine store sits inside a dense cluster of nearby stores and addresses, while a bad geocode sits alone, far from any cluster — it is a spatial outlier, and outliers are detectable without a human ever opening a map.

Prerequisites #

• Python packages: geopandas, scikit-learn (for DBSCAN and NearestNeighbors), numpy, and pandas. Install with pip install geopandas scikit-learn numpy pandas.
• A table of geocoded points with a CRS. This task assumes coordinates have already passed the structural checks — valid ranges, no nulls, parseable — described in automating coordinate validation with Python and Shapely. Spatial-outlier detection is the layer after range validation: it catches coordinates that are individually plausible but collectively wrong.
• A projected, equal-area CRS. DBSCAN’s eps and every nearest-neighbour distance must be in metres, which means projecting out of EPSG:4326 first. Use EPSG:5070 (North America Albers) for US retail data. This is non-negotiable — running DBSCAN with eps in degrees is the dominant failure mode (see below).
• An expectation of where points should be — a region, state, or market boundary — so distance-to-expected-region can flag a point that landed in the wrong part of the country entirely.

Detection Parameters #

Three complementary signals catch different failure modes; tuning their thresholds is the whole task.

eps is the radius that defines “dense enough to be normal.” For a national retail chain, 25 km is a reasonable neighbourhood scale; for a dense metro dataset, drop it. min_samples controls how many points must fall within eps for a dense core to form. Any point DBSCAN cannot assign to a dense group gets label -1 — noise — and those are the first-pass suspects.

The Nearest-Neighbour Outlier Score #

DBSCAN gives a binary in-cluster / noise label. For a graded score, measure each point’s distance to its k-th nearest neighbour and express how far that sits from the dataset’s typical spacing. For point $i$ let $d_{i,k}$ be the distance (in metres) to its $k$-th nearest neighbour. With mean $\mu$ and standard deviation $\sigma$ of those distances across all points, the outlier score is the standard z-score:

$$z_i = \frac{d_{i,k} - \mu}{\sigma}$$

A point with $z_i > 3$ sits more than three standard deviations beyond the typical neighbour distance — its nearest real neighbours are far away, which is exactly the signature of a transposed or fallback geocode dropped into empty space. Because the distribution of neighbour distances is right-skewed, applying the z-score to $\log d_{i,k}$ is often more robust; the code below computes it on the log distance.

Choosing $k$ matters. With $k = 1$, the score reacts to a single nearest neighbour, so two bad geocodes that happen to land near each other mask one another — each looks like it has a close neighbour. Setting $k$ to 5 forces a point to have several nearby real neighbours before it counts as normal, which defeats small clusters of duplicated errors. The cost is sensitivity in genuinely sparse areas, where even correct points have distant fifth-nearest neighbours; that is precisely why the region check exists as a second, independent signal.

Annotated Implementation #

The script projects to an equal-area CRS, runs DBSCAN to label noise, computes the k-NN z-score, checks distance to each point’s expected region, and combines the three into a single flag column with a reason. Every geometric step asserts its CRS.

import numpy as np
import pandas as pd
import geopandas as gpd
from pyproj import CRS
from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors

EQUAL_AREA_CRS = CRS.from_epsg(5070)   # North America Albers; metres
STORAGE_CRS = CRS.from_epsg(4326)      # interchange CRS of the input


def detect_outliers(
    gdf: gpd.GeoDataFrame,
    eps_m: float = 25_000,
    min_samples: int = 5,
    k: int = 5,
    z_threshold: float = 3.0,
    region_max_m: float = 50_000,
) -> gpd.GeoDataFrame:
    """Flag geocoding outliers using DBSCAN noise labels, a k-NN z-score, and
    distance to each point's expected region. Returns the frame with added
    'flag' (bool) and 'flag_reason' columns."""
    # --- CRS discipline: refuse to proceed without a known CRS, then project
    #     to an equal-area system so all distances are in metres. ---
    assert gdf.crs is not None, "input has no CRS; cannot compute metric distance"
    if CRS.from_user_input(gdf.crs) != STORAGE_CRS:
        gdf = gdf.to_crs(STORAGE_CRS)
    proj = gdf.to_crs(EQUAL_AREA_CRS)
    assert CRS.from_user_input(proj.crs) == EQUAL_AREA_CRS

    coords = np.column_stack([proj.geometry.x.values, proj.geometry.y.values])

    # --- Signal 1: DBSCAN noise label. eps is in METRES (safe: projected). ---
    labels = DBSCAN(eps=eps_m, min_samples=min_samples).fit_predict(coords)
    is_noise = labels == -1

    # --- Signal 2: k-NN distance z-score (on log distance, right-skew-robust). ---
    nn = NearestNeighbors(n_neighbors=k + 1).fit(coords)  # +1: self is neighbour 0
    dists, _ = nn.kneighbors(coords)
    d_k = dists[:, k]                                      # distance to k-th neighbour
    log_d = np.log1p(d_k)
    z = (log_d - log_d.mean()) / (log_d.std() + 1e-9)
    z_outlier = z > z_threshold

    # --- Signal 3: distance to the point's expected region centroid. ---
    #     'region_centroid' geometry must already be in the same projected CRS.
    region_dist = proj.geometry.distance(proj["region_centroid"])
    far_from_region = region_dist > region_max_m

    # --- Combine into one flag with a human-readable reason. ---
    reasons = []
    for noise, zo, far in zip(is_noise, z_outlier, far_from_region):
        parts = []
        if noise:
            parts.append("dbscan_noise")
        if zo:
            parts.append("knn_zscore")
        if far:
            parts.append("outside_region")
        reasons.append(",".join(parts))

    out = gdf.copy()
    out["knn_z"] = z
    out["flag"] = is_noise | z_outlier | far_from_region
    out["flag_reason"] = reasons
    return out


if __name__ == "__main__":
    stores = gpd.read_file("stores_geocoded.geojson")   # EPSG:4326
    flagged = detect_outliers(stores)
    suspects = flagged[flagged["flag"]]
    print(f"flagged {len(suspects)} of {len(flagged)} points")
    print(suspects[["store_id", "knn_z", "flag_reason"]].head(20))

The three signals are deliberately redundant. A transposed lat/lon usually trips all three — it lands in empty space (DBSCAN noise), far from any neighbour (high z), and outside its region. A ZIP-centroid fallback may only trip the k-NN score, because it is close to the right area but subtly offset from the real cluster of rooftop geocodes. Combining them catches more than any single test and, via flag_reason, tells a reviewer why each point was flagged.

The distinct signatures are worth internalizing, because they map to distinct geocoder faults. A hemisphere error — a dropped sign on longitude that puts a US store in western China, or a dropped sign on latitude that flips it into the southern hemisphere — is the most extreme case: it trips all three signals with a huge z-score and a region distance of thousands of kilometres. A lat/lon transposition typically lands the point in a different but still plausible-looking part of the map, so it reliably trips the region check and usually the k-NN score. A ZIP or city centroid fallback is the subtlest: the point is in the right metro, so the region check passes and DBSCAN may still assign it to a dense group, but its distance to the true rooftop location nudges the k-NN z-score up. Treating the three signals as a menu rather than a single test is what lets one detector span this whole range of severities.

Failure Modes and Debugging #

The dominant failure is CRS in degrees breaking eps. DBSCAN treats coordinates as points in a plane and eps as a plane distance. Feed it EPSG:4326 longitude/latitude and one “unit” of eps is one degree — roughly 111 km at the equator but varying with latitude — so the clustering is both wrong and inconsistent across the country. Projecting to an equal-area metric CRS first makes eps = 25_000 mean exactly 25 km everywhere. The second common trap is sparse regions: a legitimately remote store has no dense neighbourhood and DBSCAN labels it noise. The k-NN z-score and the region check are the counterweight — a remote-but-correct store is far from neighbours yet still inside its expected region, so requiring corroborating signals before hard-rejecting a point avoids discarding good rural sites.

Verification #

Confirm the detector is calibrated, not just running:

1. Flagged count is plausible. A healthy batch flags a small single-digit percentage. Flagging 40% means eps or the threshold is wrong (likely the degrees-vs-metres trap); flagging 0% means the threshold is too loose.
2. Manual spot-check by bounding box. Take the flagged points and inspect their bounding box against the expected footprint — flagged points should sit in empty space or the wrong region, not scatter through dense areas.
3. Reasons make sense. Group by flag_reason: transposition-style errors should show all three reasons; ZIP-centroid fallbacks typically only knn_zscore.
4. Known-bad injection. Insert a deliberately transposed copy of a real store (swap x and y) and confirm the detector flags it. A test point it misses means the thresholds are too permissive.
5. CRS asserted end to end. The two assertions in detect_outliers must pass; a projected frame that is still in degrees is the failure this guards against.

The bounding-box spot-check is a two-line inspection:

suspects = flagged[flagged["flag"]]
print("flagged bbox:", suspects.total_bounds)   # minx, miny, maxx, maxy
# A flagged point inside the main dense bbox warrants re-tuning.

Once tuned, this check becomes a standing gate rather than a one-off cleanup. Wiring it into a scheduled pipeline lets it catch geocoder regressions as new data arrives — the same principle as monitoring spatial data drift in CI pipelines, where a sudden jump in the flagged-outlier rate signals that an upstream geocoding provider changed behaviour.

Related #

• Data Validation Rules for Store Coordinates — the broader coordinate-quality regime this outlier check plugs into.
• Automating Coordinate Validation with Python and Shapely — the structural range checks that run before density-based detection.
• Monitoring Spatial Data Drift in CI Pipelines — turning a one-off check into a continuous quality gate.

← Back to Data Validation Rules for Store Coordinates