Handling Spatial Autocorrelation

Detect and handle spatial autocorrelation in ML: quantify spatial dependence with Moran's I, engineer topology-aware features, prevent data leakage, and monitor spatial drift in production.

Spatial autocorrelation — the statistical tendency for observations in close geographic proximity to exhibit similar attribute values — fundamentally violates the independent and identically distributed (IID) assumption underlying most machine learning algorithms. When training predictive models on geospatial data, ignoring this property inflates validation metrics, masks data leakage, and produces brittle inference pipelines that fail when deployed to new regions.

The production-ready response is a four-stage workflow: quantify spatial dependence, engineer topology-aware features, enforce geographic separation during validation, and monitor spatial drift after deployment. This guide implements all four stages end-to-end. It is part of Training Geospatial Predictive Models in Python, which covers the full lifecycle from feature engineering through inference automation. For a deep dive on the validation side of this workflow, see Spatial Cross-Validation Strategies.

Spatial Autocorrelation Handling Workflow Four sequential pipeline stages: Quantify Dependence (Moran's I), Engineer Features (spatial lag, ESF, decay), Spatially Blocked Validation, and Drift Monitoring. Each stage feeds into the next via arrows. 1. Quantify Dependence Moran's I · LISA 2. Engineer Topology Features lag · ESF · decay 3. Spatial Block Validation GroupKFold · buffer 4. Drift Monitoring residual I · retrain 1 2 3 4 Production pipeline for spatial dependence control

Prerequisites & Environment Setup

Before implementing spatial autocorrelation controls, your environment and data must meet these baseline requirements. Pipelines fail silently when CRS assumptions or library versions are violated.

# Pinned requirements for spatial autocorrelation handling
geopandas==0.14.4
libpysal==4.10.0
esda==2.6.0
scikit-learn==1.5.1
numpy==1.26.4
shapely==2.0.5
pyproj==3.6.1
xarray==2024.3.0
rioxarray==0.15.5
scipy==1.13.1

Install with:

pip install "geopandas==0.14.4" "libpysal==4.10.0" "esda==2.6.0" \
            "scikit-learn==1.5.1" "numpy==1.26.4" "shapely==2.0.5" \
            "pyproj==3.6.1" "xarray==2024.3.0" "rioxarray==0.15.5" "scipy==1.13.1"

GDAL/PROJ system dependencies (Ubuntu/Debian):

sudo apt-get install -y gdal-bin libgdal-dev libproj-dev

Hard prerequisites

  • Projected CRS — distance-based weight matrices require linear units (metres or feet). Geographic coordinates (latitude/longitude) produce distorted distances at higher latitudes and must be reprojected before any spatial statistics step. Full reprojection guidance is in CRS Alignment and Projection Handling.
  • Valid vector topology — no self-intersections, unclosed rings, or degenerate geometries. Run gdf.geometry.is_valid.all() and fix failures with shapely.make_valid before constructing weights.
  • Sufficient memory for sparse matrices — spatial weight matrices scale sub-quadratically when using sparse representations, but dense matrices will exhaust RAM on datasets above ~50 k observations. Always use libpysal’s sparse storage.

Validate both prerequisites at ingestion:

import geopandas as gpd

gdf = gpd.read_file("training_data.gpkg")
assert gdf.crs is not None, "GeoDataFrame has no CRS"
assert gdf.crs.is_projected, f"Expected projected CRS, got {gdf.crs.to_string()}"
assert gdf.geometry.is_valid.all(), "Invalid geometries detected — run shapely.make_valid"

Step-by-Step Implementation

Step 1 — Quantify spatial dependence with Global Moran’s I

The first step in any geospatial ML pipeline is measuring the degree of spatial clustering. Global Moran’s I remains the industry standard for detecting overall spatial autocorrelation, while Local Indicators of Spatial Association (LISA) identify hotspots and cold spots that may bias model training.

The esda library, part of the PySAL ecosystem, provides both. Always use permutation-based inference (≥ 999 permutations) rather than analytical approximations — real-world spatial boundaries rarely meet the distributional assumptions of the theoretical test.

import geopandas as gpd
import libpysal
from esda.moran import Moran, Moran_Local
import numpy as np

# Load and project to metric CRS
gdf = gpd.read_file("training_data.gpkg")
if gdf.crs.is_geographic:
    gdf = gdf.to_crs("EPSG:3857")

# Construct Queen contiguity weights (polygons)
# Use from_geodataframe for modern libpysal (>= 4.6) compatibility
w = libpysal.weights.Queen.from_geodataframe(gdf)
w.transform = "r"   # row-standardise so each row sums to 1

# ── Global Moran's I ──────────────────────────────────────────────────────
target_var = "yield_per_hectare"
y = gdf[target_var].values

moran = Moran(y, w, permutations=9999)
print(f"Global Moran's I : {moran.I:.4f}")
print(f"Expected (random): {moran.EI:.4f}")
print(f"p-value (sim)    : {moran.p_sim:.4f}")

assert moran.p_sim is not None, "Permutation test did not run — check w and y shapes"

A significant positive Moran’s I (p < 0.05) confirms spatial dependence. Values near +1 indicate strong clustering; values near 0 suggest spatial randomness; negative values imply dispersion. For raster workflows, convert pixel centroids to point geometries and use libpysal.weights.DistanceBand or KNN weights.

The Spatial Lag and Neighborhood Statistics section covers the underlying spatial weights theory, and Computing Local Moran’s I for Feature Engineering provides a full LISA implementation for identifying hotspot and coldspot geometry that would otherwise bias training.

Step 2 — Engineer topology-aware features

Once dependence is quantified, you must model it explicitly rather than hoping the algorithm absorbs it as noise. Three techniques cover the main use cases, and they are complementary rather than mutually exclusive.

2a — Spatial lag features

A spatial lag is the weighted average of neighbouring values: W * X. Appending it to your feature matrix gives gradient boosting and tree ensemble models explicit context about local neighbourhood states. It is particularly effective for environmental covariates like soil moisture, elevation gradients, and urban density where local averaging aligns with physical processes.

import numpy as np
import libpysal
import geopandas as gpd

def add_spatial_lag_features(
    gdf: gpd.GeoDataFrame,
    feature_cols: list[str],
    w: libpysal.weights.W,
) -> gpd.GeoDataFrame:
    """Append W*X spatial lag columns to a GeoDataFrame.

    Args:
        gdf: GeoDataFrame in a projected metric CRS.
        feature_cols: Names of numeric columns to lag.
        w: Row-standardised spatial weights matrix.
    Returns:
        GeoDataFrame with new '<col>_lag' columns appended.
    """
    gdf = gdf.copy()
    for col in feature_cols:
        lag_values = libpysal.weights.lag_spatial(w, gdf[col].values)
        gdf[f"{col}_lag"] = lag_values
    return gdf

feature_cols = ["elevation", "ndvi", "soil_moisture"]
gdf_lagged = add_spatial_lag_features(gdf, feature_cols, w)
assert all(f"{c}_lag" in gdf_lagged.columns for c in feature_cols), \
    "Lag columns not created — check feature_cols list"
print(f"Added {len(feature_cols)} spatial lag features")

2b — Eigenvector spatial filtering (ESF)

ESF decomposes the weights matrix into orthogonal eigenvectors representing distinct spatial scales of autocorrelation. Including the top eigenvectors as covariates effectively “whitens” residuals, restoring the IID assumption for linear models and interpretability-focused workflows. This is the preferred approach when you need statistically valid coefficient estimates rather than raw predictive performance.

from scipy.linalg import eigh
import numpy as np

def compute_spatial_eigenvectors(
    w: libpysal.weights.W,
    n_vectors: int = 10,
) -> np.ndarray:
    """Extract the top eigenvectors of the spatial weights matrix for ESF.

    Args:
        w: Row-standardised spatial weights (libpysal W object).
        n_vectors: Number of eigenvectors to retain (spatial scales to model).
    Returns:
        Array of shape (n_obs, n_vectors) — append as columns to X.
    """
    w_dense = w.full()[0]
    # Symmetrise: (W + W.T) / 2 ensures real eigenvalues
    w_sym = (w_dense + w_dense.T) / 2.0
    n = w_sym.shape[0]
    # Compute top n_vectors eigenvectors (largest eigenvalues = broadest spatial patterns)
    eigenvalues, eigenvectors = eigh(w_sym, subset_by_index=[n - n_vectors, n - 1])
    # Sort descending
    order = np.argsort(-eigenvalues)
    return eigenvectors[:, order]

esf_vecs = compute_spatial_eigenvectors(w, n_vectors=10)
assert esf_vecs.shape == (len(gdf), 10), "ESF output shape mismatch"
print(f"Extracted {esf_vecs.shape[1]} spatial eigenvectors")

2c — Distance decay features

Binary contiguity weights treat all adjacent polygons equally. Replacing them with inverse-distance or exponential decay kernels better reflects Tobler’s First Law of Geography and reduces edge artefacts in irregular tessellations.

import libpysal
import numpy as np

def build_distance_decay_weights(
    gdf: gpd.GeoDataFrame,
    bandwidth: float,
    decay: str = "exponential",
) -> libpysal.weights.W:
    """Build distance-decay spatial weights with exponential or inverse kernel.

    Args:
        gdf: GeoDataFrame in a projected metric CRS.
        bandwidth: Decay half-distance in CRS units (metres for UTM).
        decay: 'exponential' (exp(-d/bandwidth)) or 'inverse' (1/(1 + d/bandwidth)).
    Returns:
        Row-standardised libpysal W object.
    """
    coords = list(zip(gdf.geometry.x, gdf.geometry.y))
    # DistanceBand captures all neighbours within 3× bandwidth; prune with kernel
    w_raw = libpysal.weights.DistanceBand.from_array(
        np.array(coords), threshold=bandwidth * 3.0, binary=False
    )
    if decay == "exponential":
        for i in w_raw.neighbors:
            for j_idx, j in enumerate(w_raw.neighbors[i]):
                d = np.linalg.norm(
                    np.array(coords[i]) - np.array(coords[j])
                )
                w_raw.weights[i][j_idx] = np.exp(-d / bandwidth)
    w_raw.transform = "r"
    return w_raw

When integrating these engineered features into standard ML workflows, ensure your transformers preserve spatial indices. The Training with Scikit-Learn-Geo section provides pipeline templates that wrap spatial feature extraction with scikit-learn’s ColumnTransformer and Pipeline abstractions.

Step 3 — Validate with spatially blocked folds

Random k-fold cross-validation assumes data points are exchangeable. In geospatial contexts, this assumption guarantees optimistic bias: training and validation folds share spatial neighbourhoods, so autocorrelated signal bleeds across fold boundaries. The detailed mechanics of constructing these blocked splits — including variogram-range estimation and balanced fold assignment — are covered in Spatial Cross-Validation Strategies. Here the focus is on wiring the split correctly after spatial autocorrelation controls have been applied.

The key integration point is that your spatial lag and ESF columns must be computed on the full dataset before splitting, because lag values depend on global neighbour relationships. Recomputing them inside each fold would use validation-set neighbours during training, reintroducing the leakage you are trying to prevent.

import numpy as np
import geopandas as gpd
from sklearn.model_selection import GroupKFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor

# Assume gdf_lagged has spatial lag columns and a 'block_id' column
# from the spatial blocking step (see Spatial Cross-Validation Strategies)
feature_cols_full = (
    ["elevation", "ndvi", "soil_moisture"]          # original features
    + ["elevation_lag", "ndvi_lag", "soil_moisture_lag"]  # spatial lags (Step 2a)
)

X = gdf_lagged[feature_cols_full].values
# Append ESF eigenvectors
X = np.hstack([X, esf_vecs])
y = gdf_lagged["yield_per_hectare"].values
groups = gdf_lagged["block_id"].values   # geographic block IDs

estimator = GradientBoostingRegressor(
    n_estimators=300, learning_rate=0.05, max_depth=4, random_state=42
)
gkf = GroupKFold(n_splits=5)

scores = cross_val_score(
    estimator, X, y,
    cv=gkf,
    groups=groups,
    scoring="neg_root_mean_squared_error",
    n_jobs=-1,
)
print(f"Spatial CV RMSE : {-scores.mean():.3f} ± {scores.std():.3f}")
print(f"Per-fold RMSE   : {(-scores).round(3).tolist()}")

Also see Reducing Spatial Leakage in Model Training for targeted fixes covering the most common residual leakage patterns in this step.

Step 4 — Monitor spatial drift in production

Spatial autocorrelation is not static. It evolves with land-use changes, climate shifts, and infrastructure development. A model trained in 2022 may degrade by 2025 not because of algorithmic failure but because the underlying spatial dependence structure has shifted. Monitoring must track the structure of spatial signal, not just aggregate accuracy.

import numpy as np
import libpysal
from esda.moran import Moran
import geopandas as gpd
from datetime import datetime

def compute_spatial_drift_metrics(
    gdf_inference: gpd.GeoDataFrame,
    predictions: np.ndarray,
    actuals: np.ndarray,
    w_training: libpysal.weights.W,
    baseline_I: float,
    drift_threshold: float = 0.15,
) -> dict:
    """Compute spatial drift indicators for a deployed model.

    Compares current Moran's I on residuals against a training-time baseline.
    Returns a dict with drift flag, current I, delta, and a retraining recommendation.

    Args:
        gdf_inference: GeoDataFrame of inference observations (projected CRS).
        predictions: Model predictions aligned with gdf_inference rows.
        actuals: Ground-truth values for the inference batch.
        w_training: Spatial weights used during training (row-standardised).
        baseline_I: Moran's I of residuals measured at training time.
        drift_threshold: Absolute change in I that triggers a retraining alert.
    Returns:
        Dict with keys: timestamp, moran_I, delta_I, drift_detected, recommend_retrain.
    """
    residuals = actuals - predictions
    # Reconstruct weights for inference geometry (may differ from training extent)
    w_inf = libpysal.weights.Queen.from_geodataframe(gdf_inference)
    w_inf.transform = "r"
    moran_resid = Moran(residuals, w_inf, permutations=999)

    delta = abs(moran_resid.I - baseline_I)
    return {
        "timestamp": datetime.utcnow().isoformat(),
        "moran_I": round(moran_resid.I, 4),
        "p_sim": round(moran_resid.p_sim, 4),
        "delta_I": round(delta, 4),
        "drift_detected": delta > drift_threshold,
        "recommend_retrain": delta > drift_threshold and moran_resid.p_sim < 0.05,
    }

Inference pipelines should also:

  1. Validate incoming geometries against the training CRS and topology rules.
  2. Reconstruct spatial weights dynamically using cached or streaming neighbour indices.
  3. Apply spatial lag transformations consistently across batch and real-time endpoints.
  4. Flag predictions in regions with high spatial uncertainty or extrapolation risk.

When recommend_retrain is True, initiate targeted retraining with recent spatial samples rather than a full retrain — the spatial dependence shift is typically localised.

Verification & Testing

After running the full pipeline, verify each stage independently before promoting to production.

import numpy as np
from esda.moran import Moran

def verify_autocorrelation_pipeline(
    gdf,
    w,
    target_col: str,
    lag_cols: list[str],
    esf_matrix: np.ndarray,
    cv_scores: np.ndarray,
) -> None:
    """Run assertions across all four pipeline stages.

    Raises AssertionError with a descriptive message if any check fails.
    """
    # Stage 1: Moran's I ran and returned a valid statistic
    moran = Moran(gdf[target_col].values, w, permutations=999)
    assert -1.0 <= moran.I <= 1.0, f"Moran's I out of range: {moran.I}"
    print(f"[Stage 1] Moran's I = {moran.I:.4f}, p = {moran.p_sim:.4f}")

    # Stage 2: spatial lag columns are present and non-constant
    for col in lag_cols:
        assert col in gdf.columns, f"Lag column '{col}' missing from GeoDataFrame"
        assert gdf[col].std() > 0, f"Lag column '{col}' is constant — weights may be all-zero"
    print(f"[Stage 2] {len(lag_cols)} lag columns verified")

    # Stage 2b: ESF matrix shape and orthogonality check
    n, k = esf_matrix.shape
    assert n == len(gdf), "ESF row count does not match GeoDataFrame length"
    gram = esf_matrix.T @ esf_matrix
    off_diag = gram - np.diag(np.diag(gram))
    assert np.abs(off_diag).max() < 1e-6, "ESF eigenvectors are not orthogonal"
    print(f"[Stage 2] ESF orthogonality verified ({k} vectors)")

    # Stage 3: spatial CV RMSE is finite and positive
    assert np.all(np.isfinite(cv_scores)), "CV scores contain NaN or Inf"
    assert np.all(-cv_scores > 0), "Negative RMSE values — check scoring sign"
    print(f"[Stage 3] Spatial CV RMSE = {-cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

    print("[ALL STAGES] Pipeline verification passed")

Numerical sanity check — compare spatial CV against random k-fold. A substantial gap confirms that spatial autocorrelation was present and is now being controlled:

from sklearn.model_selection import KFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor

random_scores = cross_val_score(
    GradientBoostingRegressor(n_estimators=300, learning_rate=0.05, random_state=42),
    X, y,
    cv=KFold(n_splits=5, shuffle=True, random_state=42),
    scoring="neg_root_mean_squared_error",
)
spatial_rmse = -scores.mean()
random_rmse  = -random_scores.mean()
inflation_pct = (random_rmse - spatial_rmse) / spatial_rmse * 100

print(f"Random CV RMSE  : {random_rmse:.3f}")
print(f"Spatial CV RMSE : {spatial_rmse:.3f}")
print(f"Optimism gap    : {inflation_pct:.1f}%")
if inflation_pct > 10:
    print("  → Significant leakage in random CV; spatial controls are effective")
elif inflation_pct < 2:
    print("  → Low autocorrelation, or blocks may still be too small")

Troubleshooting & Common Errors

libpysal.weights.Queen.from_geodataframe raises ValueError: Island(s) detected

Cause: One or more polygons share no boundary with any neighbour (islands), leaving them with zero weight entries. This causes w.transform = "r" to fail for those rows.

Fix: Remove island observations before fitting weights, or use libpysal.weights.fill_diagonal to assign a self-loop with a small weight:

islands = [i for i, nbrs in w.neighbors.items() if len(nbrs) == 0]
print(f"Removing {len(islands)} island geometries")
gdf_clean = gdf.drop(index=islands).reset_index(drop=True)
w = libpysal.weights.Queen.from_geodataframe(gdf_clean)
w.transform = "r"

esda.moran.Moran returns NaN for I

Cause: The target variable has zero variance (all values identical) or contains NaN values. Moran’s I is undefined for a constant attribute.

Fix:

assert gdf[target_var].notna().all(), "NaN values in target — impute before computing Moran's I"
assert gdf[target_var].std() > 0, "Target variable is constant — Moran's I undefined"

MemoryError when building spatial weights on large datasets (> 200 k observations)

Cause: Dense weight matrix representation exceeds available RAM. Full NxN storage requires N² floats.

Fix: Use libpysal.weights.KNN with a small k (5–15 neighbours) rather than contiguity weights, and ensure sparse storage is active:

import libpysal, numpy as np

coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
w = libpysal.weights.KNN.from_array(coords, k=8)
w.transform = "r"
# Verify sparse backend
assert hasattr(w, "sparse"), "libpysal did not use sparse storage"

Spatial lag columns are identical for all observations

Cause: The weights matrix was not row-standardised (w.transform = "r" not called) and libpysal.weights.lag_spatial returns the raw sum of neighbour values. When all polygons have the same number of neighbours this produces a different scale, not identical values — but for point-in-polygon joins where every point has k=1 neighbours it creates a copy of the original column.

Fix: Always call w.transform = "r" immediately after constructing weights, before any lag_spatial call. Verify with assert w.transform == "R".

RuntimeError: Singular matrix encountered during ESF decomposition

Cause: The symmetrised weight matrix is singular or near-singular, typically because the weights are extremely sparse (many isolated observations) or the dataset is very small (< 30 observations).

Fix: Increase k in KNN weights (minimum 3–5 neighbours per observation) and check that n_vectors < n_observations. Reduce n_vectors proportionally for small datasets:

n_vectors = min(10, len(gdf) // 10)   # at most 10% of observations

AssertionError: ESF eigenvectors are not orthogonal

Cause: Floating-point accumulation in scipy.linalg.eigh for very large or poorly conditioned matrices. Rare but possible on datasets > 50 k with extreme weight asymmetry.

Fix: Symmetrise more aggressively and use driver="evd" for better numerical stability:

from scipy.linalg import eigh
w_sym = (w_dense + w_dense.T) / 2.0
np.fill_diagonal(w_sym, 0)   # zero self-loops
eigenvalues, eigenvectors = eigh(w_sym, subset_by_index=[n - n_vectors, n - 1], driver="evd")

Performance Optimisation

Use sparse KNN weights for large datasets

Replace contiguity weights with KNN (k=8–15) for datasets above 50 k observations. Sparse KNN weights reduce memory from O(N²) to O(N·k) and cut lag_spatial computation from O(N²) to O(N·k):

import libpysal, numpy as np

def build_knn_weights(gdf: gpd.GeoDataFrame, k: int = 10) -> libpysal.weights.W:
    coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
    w = libpysal.weights.KNN.from_array(coords, k=k)
    w.transform = "r"
    return w

Parallelise ESF with partial eigendecomposition

For large N, computing all N eigenvectors is unnecessary. Use scipy.sparse.linalg.eigsh to compute only the top k eigenvectors in O(N·k) time rather than O(N³):

from scipy.sparse.linalg import eigsh
from scipy.sparse import csr_matrix

w_sparse = csr_matrix(w.full()[0])
w_sym_sparse = (w_sparse + w_sparse.T) / 2.0
eigenvalues, eigenvectors = eigsh(w_sym_sparse, k=n_vectors, which="LM")
order = np.argsort(-eigenvalues)
esf_vecs = eigenvectors[:, order]

This reduces ESF computation from minutes to seconds on 100 k-row datasets.

Cache spatial weights between runs

Spatial weights are deterministic for a given geometry set. Persist them to avoid recomputing at each pipeline invocation:

import hashlib, pickle
from pathlib import Path

def cached_weights(gdf: gpd.GeoDataFrame, k: int = 10, cache_dir: str = ".cache"):
    key = hashlib.md5(
        str(sorted([(g.x, g.y) for g in gdf.geometry])).encode()
    ).hexdigest()[:12]
    path = Path(cache_dir) / f"weights_{key}_k{k}.pkl"
    if path.exists():
        with open(path, "rb") as f:
            return pickle.load(f)
    w = build_knn_weights(gdf, k=k)
    path.parent.mkdir(exist_ok=True)
    with open(path, "wb") as f:
        pickle.dump(w, f)
    return w

Compute spatial lag in batch chunks with Dask

For continental-scale datasets, split the GeoDataFrame into spatial partitions and compute lag features in parallel using dask-geopandas:

import dask_geopandas as dgpd
import numpy as np

dgdf = dgpd.from_geopandas(gdf, npartitions=8)
# Partition by spatial tiles so neighbours stay within partitions
dgdf_tiles = dgdf.spatial_shuffle(by="hilbert")
# Build intra-partition weights and lag per partition
# (cross-partition neighbours introduce small boundary error — acceptable for large N)

For memory-efficient distance matrix construction that feeds directly into lag computations, see Vector Proximity and Buffer Generation.

FAQ

How do I know if spatial autocorrelation is hurting my model?

Compare RMSE under random k-fold against spatially blocked cross-validation using the same estimator. If random CV reports substantially lower error — a common finding is 15–40 % lower RMSE — the gap represents leakage from autocorrelated neighbours in the training fold. Also compute Moran’s I on model residuals after fitting: a significant positive value (I > 0.1, p < 0.05) means the model has not absorbed the spatial structure and there is structured signal left in the errors.

What is the difference between spatial lag features and eigenvector spatial filtering?

Spatial lag features (W * X) add a neighbourhood average as a direct covariate. They are fast, interpretable, and work well with tree ensembles and gradient boosting. Eigenvector spatial filtering (ESF) decomposes the weights matrix into orthogonal basis functions representing spatial scales, then adds them as covariates to whiten residuals. ESF is preferred for linear models, generalised linear models, and any setting where you need statistically valid standard errors — the eigenvectors absorb spatial structure without adding correlated columns the way lag features do.

Does spatial autocorrelation affect classification as well as regression?

Yes. In classification, neighbouring observations tend to share class labels (land cover patches, urban zones, soil types). Random splits place correlated observations in both training and validation, allowing the model to memorise spatial patches rather than learning decision boundaries that generalise across regions. Apply the same spatially blocked validation regardless of task type. The only difference is that Moran’s I is computed on the predicted class probabilities or a numeric encoding of the target, not on continuous residuals.

How often should I recheck spatial autocorrelation in production?

Recompute Moran’s I on your target variable and model residuals whenever you ingest a new batch of training data, and on a scheduled cadence (monthly or quarterly) for deployed inference pipelines. Significant shifts in Moran’s I or in spatial weight matrix density signal that the underlying spatial structure has changed — this commonly happens after land-use transitions, infrastructure construction, or climate anomalies — and the model’s topology-aware features may be misaligned with the new spatial regime.

Part of: Training Geospatial Predictive Models in Python