Spatial Autocorrelation Metrics: Global & Local Statistics in Python

Spatial autocorrelation metrics quantify how strongly the values of a geographic variable resemble — or diverge from — those at neighbouring locations. They are the primary diagnostic for detecting non-random spatial structure in areal and point-referenced data, rooted in Tobler’s First Law of Geography. When ignored, spatial dependence violates the independence assumption of classical statistics, inflating Type I error rates and degrading predictive models. This page covers the full workflow: theory, construction of spatial weight matrices, global and local statistics, and production-grade validation — all within the context of the broader Core Concepts of Spatial Statistics & Geostatistics framework.

Prerequisites #

• Python 3.10+
• geopandas>=1.0, libpysal>=4.8, esda>=2.5
• numpy>=1.24, scipy>=1.10, matplotlib>=3.7, statsmodels>=0.14
• A GeoDataFrame with a single numeric attribute column and a projected CRS (e.g., UTM or LAEA — not EPSG:4326)
• Familiarity with contiguity-based and distance-based weight definitions (see spatial weight matrices)

Mathematical Core #

Global Moran’s I #

$$I = \frac{n}{\sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij}(x_i - \bar{x})(x_j - \bar{x})}{\sum_i (x_i - \bar{x})^2}$$

Where $n$ is the number of spatial units, $w_{ij}$ is the weight between units $i$ and $j$ (from the row-standardised matrix $\mathbf{W}$), $x_i$ is the attribute value at location $i$, and $\bar{x}$ is the global mean. The expected value under spatial randomness is:

$$E[I] = -\frac{1}{n-1}$$

Values significantly above $E[I]$ indicate positive autocorrelation (clustering); values below indicate negative autocorrelation (dispersion).

The expression is easier to trust once you read it as two factors. The right-hand ratio is a covariance-to-variance comparison: the numerator accumulates the product of paired deviations $(x_i - \bar{x})(x_j - \bar{x})$ over every neighbour pair, so it is positive whenever two neighbours sit on the same side of the mean and negative whenever they straddle it, while the denominator is the total variance, which strips the units out of the result. The left-hand factor $n / \sum_i \sum_j w_{ij}$ rescales for the size of the weight matrix. Under row standardisation every row sums to 1, so $\sum_i \sum_j w_{ij} = n$ and that factor collapses to 1 — which is precisely why row-standardised weights let $I$ be read as the slope of spatial lag regressed on value. Note also that $I$ is not bounded to $[-1, 1]$: it is a ratio of quadratic forms, and its attainable extremes are set by the smallest and largest eigenvalues of the symmetrised weight matrix. On irregular polygon layouts with very uneven neighbour counts, a legitimate $I$ can exceed 1 in magnitude.

Geary’s C #

$$C = \frac{(n-1)}{2\sum_i \sum_j w_{ij}} \cdot \frac{\sum_i \sum_j w_{ij}(x_i - x_j)^2}{\sum_i (x_i - \bar{x})^2}$$

Geary’s C operates on pairwise squared differences rather than deviations from the mean, making it more sensitive to local divergences. It ranges from 0 to 2: values below 1 indicate positive autocorrelation; values above 1 indicate negative autocorrelation.

The numerator $\sum_i \sum_j w_{ij}(x_i - x_j)^2$ is the same quantity that drives the empirical semivariogram, so Geary’s C behaves like a variogram evaluated over a single neighbour lag rather than like a correlation. That relationship explains the difference in sensitivity: one extreme value adjacent to a typical one contributes a large squared difference to $C$ but only a moderate cross-product to $I$, so $C$ moves further for the same anomaly. The two statistics are related without being redundant — the informal rule $C \approx 1 - I$ holds only when the variable is roughly symmetric and neighbour counts are close to uniform. A wide gap between them is itself diagnostic, and usually points to skew or to a handful of dominant outliers rather than to a disagreement about whether structure exists.

Local Moran’s I (LISA) #

$$I_i = \frac{(x_i - \bar{x})}{m_2} \sum_j w_{ij}(x_j - \bar{x}), \quad m_2 = \frac{\sum_i (x_i - \bar{x})^2}{n}$$

Each feature receives its own statistic $I_i$, enabling classification into high-high, low-low, high-low, and low-high spatial association quadrants.

Workflow Diagram #

Step 1: Data Preparation & Weight Construction #

Spatial autocorrelation metrics are entirely dependent on the structure of the underlying weight matrix. Poorly constructed weights produce misleading coefficients regardless of which test follows.

import geopandas as gpd
import libpysal

# 1. Load and validate CRS — weights fail on unprojected lat/lon
gdf = gpd.read_file("data/urban_indicators.shp")
if not gdf.crs.is_projected:
    gdf = gdf.to_crs(epsg=32633)  # UTM Zone 33N; choose zone matching your study area

# 2. Repair topology and remove features with null attributes
gdf["geometry"] = gdf.geometry.make_valid()
gdf = gdf.dropna(subset=["geometry", "income_index"]).reset_index(drop=True)

# 3. Queen contiguity — shares edges or vertices; robust for irregular polygons
w = libpysal.weights.Queen.from_dataframe(gdf)

# 4. Warn about disconnected components (islands artificially inflate I)
if w.n_components > 1:
    print(f"Warning: {w.n_components} disconnected components detected. "
          "Consider KNN fallback for isolated features.")
    # w = libpysal.weights.KNN.from_dataframe(gdf, k=4)

# 5. Always row-standardise: spatial lag becomes a weighted average,
#    keeping coefficients comparable across regions with varying neighbour counts
w.transform = "R"
print(f"Features: {w.n}, Max neighbours: {max(w.cardinalities.values())}")

Always call make_valid() before weight construction when your data originates from administrative boundary files — clipped or dissolved geometries frequently contain self-intersections that cause silent failures in Queen.from_dataframe.

Step 2: Global Spatial Autocorrelation #

Global statistics return a single coefficient describing structure across the entire study area. Use them to answer: “Is there any evidence of non-random spatial organisation?”

import esda
import numpy as np

y = gdf["income_index"].values

# --- Moran's I ---
moran = esda.Moran(y, w, permutations=999)
print(f"Moran's I:           {moran.I:.4f}")
print(f"Expected I:          {moran.EI:.4f}")        # ≈ -1/(n-1)
print(f"Permutation p-value: {moran.p_sim:.4f}")
print(f"z-score:             {moran.z_norm:.4f}")

# --- Geary's C ---
geary = esda.Geary(y, w, permutations=999)
print(f"Geary's C:           {geary.C:.4f}")         # <1 = positive autocorrelation
print(f"Geary's C p-value:   {geary.p_sim:.4f}")

Interpretation quick-reference:

Always report the permutation count alongside p-values. For publication, raise permutations to 9999 to stabilise pseudo p-values near the significance boundary.

Step 3: Moran Scatterplot #

The Moran scatterplot visualises the relationship between each feature’s standardised value ($z_i$) and its spatial lag ($\sum_j w_{ij} z_j$). The slope of the regression line equals Moran’s I. Quadrant position classifies each observation: High-High (upper-right), Low-Low (lower-left), High-Low (lower-right, spatial outlier), Low-High (upper-left, spatial outlier).

import matplotlib.pyplot as plt

z = (y - y.mean()) / y.std()
spatial_lag = libpysal.weights.lag_spatial(w, z)

fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(z, spatial_lag, alpha=0.5, s=18, color="steelblue", label="Features")
ax.axhline(0, color="gray", lw=0.8, ls="--")
ax.axvline(0, color="gray", lw=0.8, ls="--")

# Regression line (slope = Moran's I)
m = np.polyfit(z, spatial_lag, 1)
xs = np.linspace(z.min(), z.max(), 200)
ax.plot(xs, np.polyval(m, xs), color="crimson", lw=1.5,
        label=f"Moran's I = {moran.I:.3f}")

ax.set_xlabel("Standardised value (z)")
ax.set_ylabel("Spatial lag")
ax.set_title("Moran Scatterplot")
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig("moran_scatterplot.png", dpi=150)

Step 4: Local Indicators of Spatial Association (LISA) #

Global coefficients mask local heterogeneity. LISA decomposes the global signal into per-feature statistics, enabling hotspot detection and spatial outlier identification. For a step-by-step walkthrough of the Moran’s I calculation itself, see How to Calculate Moran’s I in PySAL.

from statsmodels.stats.multitest import fdrcorrection

# --- Local Moran's I ---
lisa = esda.Moran_Local(y, w, permutations=999, seed=42)

gdf["lisa_I"]     = lisa.Is        # local statistic per feature
gdf["lisa_p"]     = lisa.p_sim     # raw permutation p-value
gdf["lisa_q"]     = lisa.q         # quadrant: 1=HH, 2=LH, 3=LL, 4=HL

# Benjamini–Hochberg FDR correction — mandatory when testing hundreds of features
_, gdf["lisa_p_fdr"] = fdrcorrection(gdf["lisa_p"], alpha=0.05)
gdf["lisa_sig"]   = gdf["lisa_p_fdr"] < 0.05

# --- Getis-Ord G* ---
# Requires binary (not row-standardised) weights with self-inclusion
w_binary = libpysal.weights.Queen.from_dataframe(gdf)   # untransformed
g_star = esda.G_Local(y, w_binary, permutations=999, star=True, seed=42)

gdf["gstar_z"] = g_star.Zs        # positive z: intensity hotspot; negative: coldspot
gdf["gstar_p"] = g_star.p_sim

# Significant clusters only (FDR-adjusted)
_, gdf["gstar_p_fdr"] = fdrcorrection(gdf["gstar_p"], alpha=0.05)
gdf["gstar_sig"] = gdf["gstar_p_fdr"] < 0.05

LISA quadrant classification:

Never map raw LISA values without FDR or Bonferroni correction. Plotting significance at the raw $p < 0.05$ level across 500 features will produce roughly 25 false positives by chance alone.

LISA vs. Getis-Ord G:* LISA identifies the quadrant (direction of association), while Getis-Ord G* measures intensity (magnitude of concentration). Use them together: LISA pinpoints outliers and both-direction clusters; G* ranks the strength of hotspots and coldspots independently.

Step 5: Diagnostic Configuration #

Several validation steps must precede final interpretation.

Trend vs. True Clustering #

Before running local statistics, test whether your variable exhibits a large-scale spatial trend. Persistent trends inflate global Moran’s I and create artefact hotspots near the gradient boundaries. Apply stationarity and trend analysis — or at minimum, regress out first- and second-order polynomial surfaces — before computing LISA.

from sklearn.linear_model import LinearRegression

# Detrend using centroid coordinates
coords = np.column_stack([gdf.geometry.centroid.x, gdf.geometry.centroid.y])
X_poly = np.column_stack([coords, coords**2, coords[:,0]*coords[:,1]])
trend = LinearRegression().fit(X_poly, y)
y_detrended = y - trend.predict(X_poly)

moran_detrended = esda.Moran(y_detrended, w, permutations=999)
print(f"Moran's I (detrended): {moran_detrended.I:.4f}")

Residual Diagnostics #

If you are fitting spatial regression models (SAR, SEM, GWR), always run Moran’s I on residuals. Significant autocorrelation in residuals indicates model misspecification.

# After fitting a spatial model, test residuals
residuals = y - model.predict(X)
moran_resid = esda.Moran(residuals, w, permutations=999)
if moran_resid.p_sim < 0.05:
    print(f"Residual autocorrelation detected (I={moran_resid.I:.4f}). "
          "Consider spatial lag or spatial error specification.")

When the Exchangeability Assumption Weakens #

Permutation inference assumes the observed values are exchangeable under the null — that any value could plausibly have landed at any location. Rates computed as counts over a population break that assumption, because a tract of 200 residents produces a far noisier rate than one of 20,000, and the noisier values carry more sampling variance wherever they happen to fall. The failure is systematic rather than random: sparsely populated units generate extreme rates, extreme rates cluster wherever small denominators cluster, and Moran’s I ends up describing the geography of the denominator instead of the geography of the process. Use the Assunção–Reis empirical Bayes adjustment, which shrinks each rate toward the global rate in proportion to the size of its denominator, before testing:

# Rates: shrink toward the global rate before testing autocorrelation
moran_rate = esda.Moran_Rate(gdf["cases"].values,
                             gdf["population"].values,
                             w, permutations=999)   # adjusted=True by default
print(f"EB-adjusted Moran's I: {moran_rate.I:.4f}")

A milder version of the same degradation affects any heavily skewed variable, rate or not. A log or rank transform before testing usually moves $I$ toward the value the untransformed data would produce with a handful of extremes removed. If it does not move at all, the clustering is being driven by those extremes, and the result should be reported with that dependence stated rather than as a property of the variable as a whole.

Scale Sensitivity (MAUP) #

Autocorrelation coefficients are scale-dependent. Aggregating to coarser zones artificially inflates Moran’s I. Always document the spatial resolution and, where feasible, repeat the analysis at multiple zoning schemes to assess sensitivity.

Output Interpretation: What Good Looks Like vs. Warning Signs #

Moran scatterplot — healthy signals:

• Most points in the HH (upper-right) and LL (lower-left) quadrants with a clearly positive slope confirms meaningful clustering.
• A few HL/LH outliers are expected and theoretically informative.

LISA map — healthy signals:

• Geographically coherent patches of HH or LL significance (not scattered individual pixels).
• After FDR correction, significant count drops by less than 50% — if more than half your features lose significance, the raw pattern was largely noise.

Warning signs:

• Moran’s I > 0.8 on administrative boundary data almost always reflects MAUP artefacts, not genuine clustering.
• LISA significance concentrated exclusively at the study area perimeter indicates edge effects rather than real hotspots.
• Geary’s C near 0 with non-significant Moran’s I suggests a local pattern that cancels out globally — investigate LISA quadrants rather than reporting a null global result.

Production Considerations #

Memory and Performance #

esda.Moran_Local with 9999 permutations on 50,000 features requires holding the attribute array and the full permutation matrix in memory simultaneously. On a 16 GB machine:

# Chunked permutation for large datasets — run global first, local per chunk
import psutil
n_features = len(gdf)
# Approximate RAM requirement: n * permutations * 8 bytes (float64)
ram_gb = (n_features * 9999 * 8) / 1e9
print(f"Estimated RAM for LISA: {ram_gb:.1f} GB")

# For n > 20,000: reduce permutations or use seed-fixed reproducibility
perms = 9999 if n_features < 20_000 else 999
lisa = esda.Moran_Local(y, w, permutations=perms, seed=42)

Sparse Matrices #

libpysal returns sparse scipy.sparse weight matrices internally. For datasets above ~10,000 features, avoid converting to dense arrays:

# Efficient spatial lag computation using sparse representation
from scipy.sparse import issparse
W_sparse = w.sparse   # scipy CSR matrix
spatial_lag_fast = W_sparse.dot(y)  # O(nnz) instead of O(n²)

Parallelisation #

For repeated analysis (bootstrap, scenario comparison), parallelise across independent permutation seeds rather than within a single LISA call:

from joblib import Parallel, delayed

def run_lisa_seed(y, w, seed):
    return esda.Moran_Local(y, w, permutations=999, seed=seed).Is

results = Parallel(n_jobs=-1)(
    delayed(run_lisa_seed)(y, w, s) for s in range(50)
)

Troubleshooting #

Frequently Asked Questions #

What does a Moran’s I of 0.6 mean? A value of 0.6, significantly above the expected value ($E[I] \approx -1/(n-1) \approx 0$), indicates strong positive spatial autocorrelation: similar attribute values tend to be located near one another rather than distributed randomly across the study area.

When should I use Geary’s C instead of Moran’s I? Geary’s C operates on pairwise squared differences between neighbours, making it more sensitive to local heterogeneity and spatial outliers than Moran’s I, which is based on deviations from the global mean. Use them together: Moran’s I for a dataset-level summary, Geary’s C as a complementary check for local departures.

How many permutations should I use? At least 999 for exploratory analysis; 9999 for publication. Fewer permutations produce unstable pseudo p-values, particularly near the significance boundary (e.g., $p \approx 0.04$ can fluctuate to $p > 0.05$ with only 99 permutations).

Do I always need FDR correction for LISA? Yes. LISA runs one significance test per feature. At $\alpha = 0.05$ across 400 features, you expect 20 false positives by chance. Benjamini-Hochberg FDR correction (via statsmodels.stats.multitest.fdrcorrection) is the standard adjustment and takes one line of code.

Can Moran’s I be greater than 1? Yes. Moran’s I is a ratio of quadratic forms, not a correlation coefficient, and its attainable range is fixed by the extreme eigenvalues of the symmetrised weight matrix rather than by $[-1, 1]$. Values slightly above 1 are legitimate on irregular polygon layouts with very uneven neighbour counts. Values far above 1 almost always indicate duplicated geometries or a weight matrix that no longer matches the attribute array after a dropna or reset_index.

Next Steps #

For a focused walkthrough of computing Moran’s I step-by-step with version-pinned PySAL imports, see How to Calculate Moran’s I in PySAL. If your data consists of discrete event coordinates rather than contiguous areal units, the statistics described here will misrepresent spatial structure; Point Pattern Analysis provides the appropriate framework for that case. To understand how autocorrelation evidence feeds weight construction decisions, review spatial weight matrices and the building custom spatial weight matrices guide.

Related #

• How to Calculate Moran’s I in PySAL — step-by-step PySAL implementation with version pinning
• Spatial Weight Matrices — Queen, Rook, KNN, and distance-band construction
• Building Custom Spatial Weights Matrices — hybrid and user-defined weight strategies
• Stationarity & Trend Analysis — prerequisite detrending before local testing
• Point Pattern Analysis — alternative framework for event-based spatial data

← Back to Core Concepts of Spatial Statistics & Geostatistics