Spatial Data Leakage in Machine Learning: How to Fix It

The Lie Your Model Tells

Consider an illustrative production scenario. A team trains a churn model on 200,000 customer records, uses a random 80/20 split, and reports a test AUC of 0.91. Six months after deployment, performance in newly served regions is 0.67 — even though the preprocessing and model code are unchanged.

The hidden difference is geography. Nearby customers from the same cities, neighbourhoods, and service areas appeared in both training and test data. The evaluation therefore measured interpolation among familiar neighbours, while deployment required extrapolation to unfamiliar places.

This is spatial leakage: dependence across nearby observations makes a nominally held-out test set less independent than it appears. Research on spatially structured data shows that ordinary random cross-validation can produce optimistic performance estimates when the intended prediction domain differs geographically from the training sample. Roberts et al. review the underlying validation problem, while a later spatial k-fold study demonstrates it across real datasets.

Tobler’s First Law — and Why Your Model Does Not Know It

Waldo Tobler’s foundational 1970 paper states: "Everything is related to everything else, but near things are more related than distant things." The original publication captures the intuition behind spatial autocorrelation.

Standard supervised-learning workflows often treat rows as independent and identically distributed. Geographic observations may instead form correlated neighbourhoods, environmental zones, markets, or service regions. A random split preserves those local relationships across both sides of the evaluation boundary.

How Leakage Happens

Imagine a property model evaluated with a random split inside one city. Test properties share school catchments, transport access, flood exposure, and neighbourhood character with training properties a few streets away. The model has not seen the exact test rows, but it has seen much of their local context.

If the deployment target is a new city, that validation design answers the wrong question. The evaluation protocol must reproduce the separation between the available training geography and the intended prediction geography.

Measuring Spatial Autocorrelation with Moran’s I & LISA

Moran’s I summarizes spatial correlation by comparing how deviations of neighboring feature values covary compared to the total sample variance. Formally:

from libpysal.weights import Queen
from esda.moran import Moran, Moran_Local
import geopandas as gpd

gdf = gpd.read_file("observations.geojson")
# Generate spatial weights matrix with row-standardization
w = Queen.from_dataframe(gdf)
w.transform = "r"

# Audit model residuals for spatial clustering
residuals = y_test - model.predict(X_test)
moran = Moran(residuals, w)
print(f"Residual Moran's I: {moran.I:.4f} (p-value: {moran.p_sim:.4f})")
# If p < 0.05, the model's errors are not independent; spatial leakage has occurred.

Buffered Spatial Cross-Validation: The Dead-Zone Guarantee

Even when partitioning data into spatial blocks, points sitting right along fold boundaries still exhibit strong spatial autocorrelation. Buffered Spatial Cross-Validation enforces an exclusion buffer zone around the test fold:

import geopandas as gpd
from shapely.ops import unary_union

def buffered_spatial_train_test_split(gdf, test_polygon, buffer_distance_meters=2000):
    """Drops training observations within buffer_distance of the test boundary."""
    # Isolate test set within geographic target boundary
    test_mask = gdf.geometry.within(test_polygon)
    test_gdf = gdf[test_mask].copy()
    
    # Create exclusion dead-zone buffer around the test set
    test_boundary_union = unary_union(test_gdf.geometry)
    exclusion_zone = test_boundary_union.buffer(buffer_distance_meters)
    
    # Train set strictly excludes test points AND points within the buffer
    train_gdf = gdf[~gdf.geometry.intersects(exclusion_zone)].copy()
    
    print(f"Test samples: {len(test_gdf)} | Train samples: {len(train_gdf)} | Dropped buffer samples: {len(gdf) - len(test_gdf) - len(train_gdf)}")
    return train_gdf, test_gdf

Hierarchical H3 Hexagons & Positional Embeddings

Passing raw (latitude, longitude) into decision trees forces orthogonal axis splits that overfit bounding rectangles. Production GeoAI systems encode spatial context using discrete Uber H3 hexagons or multi-frequency sinusoidal embeddings:

import h3
import numpy as np

def add_hierarchical_spatial_features(df, lat_col="lat", lon_col="lon"):
    # Multi-resolution hexagonal spatial tokens
    df["h3_res6"] = df.apply(lambda r: h3.geo_to_h3(r[lat_col], r[lon_col], resolution=6), axis=1) # ~36 km²
    df["h3_res8"] = df.apply(lambda r: h3.geo_to_h3(r[lat_col], r[lon_col], resolution=8), axis=1) # ~0.7 km²
    
    # Spherical multi-frequency harmonic coordinate embedding (Sphere2Vec pattern)
    lat_rad = np.radians(df[lat_col].values)
    lon_rad = np.radians(df[lon_col].values)
    
    for freq in [1, 2, 4, 8, 16]:
        df[f"sin_lat_{freq}"] = np.sin(freq * lat_rad)
        df[f"cos_lat_{freq}"] = np.cos(freq * lat_rad)
        df[f"sin_lon_{freq}"] = np.sin(freq * lon_rad)
        df[f"cos_lon_{freq}"] = np.cos(freq * lon_rad)
        
    return df

Five-Step Action Plan

Frequently Asked Questions