Spatial Autocorrelation Explained

Problem statement

Tobler's first law — "near things are more related than distant things" — is usually quoted as an observation. It is better understood as a warning, because almost every statistical method you were taught assumes the opposite.

Ordinary regression assumes independent residuals. Confidence intervals assume independent observations. A t-test on 400 spatial units with strongly correlated values is not really working with 400 independent pieces of information — it might be closer to 40 — and its p-values are correspondingly too small.

Spatial autocorrelation is how you measure that dependence. Moran's I is the standard statistic, and it has a property that surprises people:

import numpy as np
import geopandas as gpd
from shapely.geometry import box
from esda.moran import Moran
from libpysal.weights import Queen, Rook

# a constructed 20 x 20 checkerboard: the most dispersed pattern there is
cells = gpd.GeoDataFrame({"value": [(r + c) % 2 for r in range(20) for c in range(20)]},
                         geometry=[box(c, r, c + 1, r + 1) for r in range(20) for c in range(20)])
checkerboard = {}
for name, build in (("Rook weights (4 edge neighbours)", Rook), ("Queen weights (8 neighbours)", Queen)):
    wc = build.from_dataframe(cells, use_index=True)
    wc.transform = "r"
    np.random.seed(0)
    checkerboard[name] = Moran(cells["value"].to_numpy(dtype="float64"), wc, permutations=999)
    print(f"  {name:<34} I = {checkerboard[name].I:+.4f}   p = {checkerboard[name].p_sim:.4f}")

  Rook weights (4 edge neighbours)   I = -1.0000   p = 0.0010
  Queen weights (8 neighbours)       I = -0.0393   p = 0.0700

Identical data. One weights definition says "the most dispersed pattern possible". The other says "indistinguishable from random". Neither is a bug.

Quick answer

Moran's I compares each value to the average of its neighbours:

from libpysal.weights import Queen
from esda.moran import Moran

w = Queen.from_dataframe(zones, use_index=True)
w.transform = "r"                                  # row-standardise

moran = Moran(zones["rate"].to_numpy(), w, permutations=999)
print(f"I = {moran.I:+.4f}  E[I] = {moran.EI:+.4f}  "
      f"z = {moran.z_sim:+.2f}  p = {moran.p_sim:.4f}")

I = +0.3455  E[I] = -0.0021  z = +14.03  p = 0.0010

E[I] is not zero. Under the null hypothesis it is −1/(n−1), which for Leeds' 488 neighbourhoods is −0.0021. Comparing I to zero rather than to E[I] is a small error at n=488 and a large one at n=20.

What Moran's I compares, weights, permutation inference

1. Understand what the statistic actually compares

Moran's I relates each unit's value to the weighted mean of its neighbours' values. Precisely, with row-standardised weights it is the regression slope of the spatial lag on the standardised value:

import numpy as np

values = zones["rate"].to_numpy()
z = (values - values.mean()) / values.std()
lag = w.sparse @ z                       # weighted mean of neighbours' z-scores

print(f"Moran's I             {moran.I:+.6f}")
print(f"OLS slope of lag on z {np.polyfit(z, lag, 1)[0]:+.6f}")
print(f"z'Wz / z'z            {(z @ lag) / (z @ z):+.6f}")
print(f"correlation(z, lag)   {np.corrcoef(z, lag)[0, 1]:+.6f}")

Moran's I             +0.345484
OLS slope of lag on z +0.345484
z'Wz / z'z            +0.345484
correlation(z, lag)   +0.563485

The first three agree exactly. The correlation does not — it is 0.563 against a slope of 0.345, because the spatial lag is an average and therefore less variable than the values it averages (sd(lag)/sd(z) = 0.613 here, and 0.563485 × 0.613 = 0.345484).

So I is the slope of the Moran scatterplot, not its correlation coefficient. The distinction is worth keeping straight: averaging always shrinks variance, so I is systematically smaller in magnitude than the correlation, and quoting one for the other overstates the pattern.

2. Recognise that the weights are half the analysis

The w object defines what "neighbour" means, and different definitions answer different questions. On the Leeds crime rates:

from libpysal.weights import KNN

for name, wi in (("Queen", Queen.from_dataframe(zones, use_index=True)),
                 ("Rook", Rook.from_dataframe(zones, use_index=True)),
                 ("KNN k=4", KNN.from_dataframe(zones, k=4)), ("KNN k=8", KNN.from_dataframe(zones, k=8))):
    wi.transform = "r"
    np.random.seed(0)
    mi = Moran(values, wi, permutations=999)
    print(f"  {name:<8}  mean nb {wi.mean_neighbors:4.1f}  I = {mi.I:+.4f}  p = {mi.p_sim:.4f}")

  Queen     mean nb  5.7  I = +0.3455  p = 0.0010
  Rook      mean nb  5.7  I = +0.3570  p = 0.0010
  KNN k=4   mean nb  4.0  I = +0.2866  p = 0.0010
  KNN k=8   mean nb  8.0  I = +0.2786  p = 0.0010

On real neighbourhoods queen and rook are nearly the same thing — boundaries almost never meet at a single point, so both give 5.7 neighbours on average — while KNN, which ignores boundaries, finds noticeably weaker autocorrelation. For a pattern whose structure operates at the scale of the neighbourhood definition, the choice decides the answer entirely, as the checkerboard from the problem statement shows:

for name, mc in checkerboard.items():              # the checkerboard from the problem statement
    print(f"  {name:<34} I = {mc.I:+.4f}  z = {mc.z_sim:+6.2f}  p = {mc.p_sim:.4f}")

  Rook weights (4 edge neighbours)   I = -1.0000  z = -26.66  p = 0.0010
  Queen weights (8 neighbours)       I = -0.0393  z =  -1.40  p = 0.0700

Rook weights see only the four edge-sharing neighbours, which on a checkerboard are all the opposite colour — perfect dispersion. Queen weights add the four diagonals, which are all the same colour, and the two effects cancel almost exactly.

Choose the weights from the process, and say which you used. A Moran's I quoted without its weights specification is not reproducible.

3. Row-standardise, and know what it does

w.transform = "r"

Row standardisation makes each unit's weights sum to 1, so the spatial lag is a weighted mean of neighbours rather than a sum. Without it, a unit with twelve neighbours contributes three times as much as one with four, purely because of how the boundaries fall.

This matters most for irregular administrative units, where neighbour counts vary enormously. For a regular grid it changes little.

4. Use the permutation p-value

moran = Moran(values, w, permutations=999)
print(f"analytical p = {moran.p_norm:.4f}   permutation p = {moran.p_sim:.4f}")

analytical p = 0.0000   permutation p = 0.0010

The analytical p-value assumes normality. The permutation ("conditional randomisation") p-value shuffles the values across the units many times and asks how often a random arrangement produces an I as extreme as the observed one. It makes no distributional assumption and is the right default.

With 999 permutations the smallest achievable p is 0.001 — so p = 0.0010 means "never matched in 999 tries", not "p equals exactly one in a thousand".

5. Interpret a non-significant result correctly

Shuffle the Leeds rates across the neighbourhoods — the same values with the geography removed — and test again:

shuffled = np.random.default_rng(1).permutation(values)
np.random.seed(0)
mr = Moran(shuffled, w, permutations=999)
print(f"  shuffled   I = {mr.I:+.4f}  E[I] = {mr.EI:+.4f}  z = {mr.z_sim:+.2f}  p = {mr.p_sim:.4f}")

  shuffled   I = +0.0116  E[I] = -0.0021  z = +0.53  p = 0.2490

That is the correct answer for data with no spatial structure: the same 488 rates, arranged at random. Note that I is not zero — it is +0.012, because any finite random arrangement has some incidental pattern. The z-score of 0.53 is what tells you it is unremarkable.

Reporting I alone, without the z-score or p-value, invites the reader to see structure in the second decimal place.

Global and local Moran's I on Leeds crime rates

The data used on this page

Every example runs on real data: the 488 Lower Layer Super Output Areas (LSOAs) of Leeds from the ONS, their Census 2021 population from Nomis, and the street-level crimes West Yorkshire Police published for June 2026 on data.police.uk. rate is crimes per 1,000 residents.

import io
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query")
features, offset = [], 0
while True:                                   # the service pages at 1,000 features
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD,LSOA21NM",
                                     "outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700).set_index("LSOA21CD").sort_index()

population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(zones.index)}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = population

west, south, east, north = zones.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                  data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if reply.status_code != 429:
                break
            time.sleep(2 ** attempt)
        reply.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "lon": float(c["location"]["longitude"]),
                    "lat": float(c["location"]["latitude"])} for c in reply.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes.lon, crimes.lat), crs=4326).to_crs(27700)
incidents = gpd.sjoin(incidents, zones[["geometry"]], predicate="within").rename(columns={"index_right": "lsoa", "LSOA21CD": "lsoa"})
zones["crimes"] = incidents["lsoa"].value_counts().reindex(zones.index, fill_value=0)
zones["rate"] = 1000 * zones["crimes"] / zones["population"]
print(f"{len(zones)} LSOAs, {zones['population'].sum():,} residents, {len(incidents):,} crimes; "
      f"rate median {zones['rate'].median():.1f}, max {zones['rate'].max():.1f} per 1,000")

488 LSOAs, 811,957 residents, 9,229 crimes; rate median 8.4, max 174.4 per 1,000

Example 1 — global Moran's I with the choices made explicit

import geopandas as gpd
import numpy as np
import pandas as pd
from esda.moran import Moran
from libpysal.weights import KNN, Queen, Rook


WEIGHT_BUILDERS = {
    "queen": lambda gdf, **kw: Queen.from_dataframe(gdf, use_index=True),
    "rook": lambda gdf, **kw: Rook.from_dataframe(gdf, use_index=True),
    "knn4": lambda gdf, **kw: KNN.from_dataframe(gdf, k=4),
    "knn8": lambda gdf, **kw: KNN.from_dataframe(gdf, k=8),
}


def global_moran(zones, column, scheme="queen", *, permutations=999, transform="r"):
    w = WEIGHT_BUILDERS[scheme](zones)
    w.transform = transform
    if w.islands:
        print(f"  warning: {len(w.islands)} unit(s) with no neighbours were excluded")

    values = zones[column].to_numpy()
    moran = Moran(values, w, permutations=permutations)

    print(f"{scheme:6} n={w.n} mean_nb={w.mean_neighbors:4.1f}  "
          f"I={moran.I:+.4f} E[I]={moran.EI:+.4f} z={moran.z_sim:+7.2f} p={moran.p_sim:.4f}")
    return moran, w


def weights_sensitivity(zones, column):
    rows = []
    for scheme in WEIGHT_BUILDERS:
        moran, w = global_moran(zones, column, scheme)
        rows.append({"scheme": scheme, "mean_nb": round(w.mean_neighbors, 1),
                     "I": round(moran.I, 4), "p": moran.p_sim})
    frame = pd.DataFrame(rows)
    print(f"\nI ranges {frame['I'].min():+.4f} to {frame['I'].max():+.4f} across weights")
    return frame


weights_sensitivity(zones, "rate")

queen  n=488 mean_nb= 5.7  I=+0.3455 E[I]=-0.0021 z= +13.60 p=0.0010
rook   n=488 mean_nb= 5.7  I=+0.3570 E[I]=-0.0021 z= +15.21 p=0.0010
knn4   n=488 mean_nb= 4.0  I=+0.2866 E[I]=-0.0021 z= +10.62 p=0.0010
knn8   n=488 mean_nb= 8.0  I=+0.2786 E[I]=-0.0021 z= +14.77 p=0.0010

I ranges +0.2786 to +0.3570 across weights

A spread of 0.078 across four definitions — contiguity against nearest neighbours — is worth reporting rather than hiding: every scheme finds clear clustering, but how much depends on what "neighbour" means. On the checkerboard the queen–rook gap alone is 0.96, which is the signal that the weights, not the data, are driving the answer.

Example 2 — the Moran scatterplot, which is where the interpretation lives

import matplotlib.pyplot as plt


def moran_scatter(zones, column, w, moran, ax=None):
    values = zones[column].to_numpy()
    z = (values - values.mean()) / values.std()
    lag = w.sparse @ z

    ax = ax or plt.subplots(figsize=(6, 6))[1]
    ax.axhline(0, color="#94a3b8", lw=1)
    ax.axvline(0, color="#94a3b8", lw=1)
    ax.scatter(z, lag, s=14, alpha=0.6, color="#0ea5e9")

    line = np.linspace(z.min(), z.max(), 10)
    ax.plot(line, moran.I * line, color="#ef4444", lw=2,
            label=f"slope = I = {moran.I:.3f}")

    for x, y, label in [(0.6, 0.9, "HH"), (-0.9, 0.9, "LH"),
                        (-0.9, -0.9, "LL"), (0.6, -0.9, "HL")]:
        ax.text(x * z.max(), y * lag.max(), label, fontsize=13,
                weight="bold", color="#64748b")

    ax.set_xlabel(f"{column} (standardised)")
    ax.set_ylabel("spatial lag of neighbours")
    ax.legend()
    return ax


moran, w = global_moran(zones, "rate")
moran_scatter(zones, "rate", w, moran)

queen  n=488 mean_nb= 5.7  I=+0.3455 E[I]=-0.0021 z= +14.45 p=0.0010

The four quadrants are the whole vocabulary of local spatial statistics:

• HH — high value, high neighbours: a hotspot
• LL — low value, low neighbours: a coldspot
• HL and LH — spatial outliers: a unit unlike its surroundings

The global I is the slope of the fitted line — one number covering four distinct kinds of behaviour, which is exactly why local statistics exist.

Example 3 — from global to local

from esda.moran import Moran_Local

def local_moran(zones, column, w, *, permutations=999, alpha=0.05):
    np.random.seed(42)                        # permutation results move between runs
    lisa = Moran_Local(zones[column].to_numpy(), w, permutations=permutations)

    labels = np.array(["not significant", "HH", "LH", "LL", "HL"])
    quadrant = np.where(lisa.p_sim <= alpha, lisa.q, 0)

    out = zones.copy()
    out["lisa_I"] = lisa.Is
    out["lisa_p"] = lisa.p_sim
    out["cluster"] = labels[quadrant]

    print(out["cluster"].value_counts().to_string())
    return out


result = local_moran(zones, "rate", w)

cluster
not significant    359
LL                  59
HH                  46
LH                  15
HL                   9

One hundred and eleven high-high units and one hundred and sixteen low-low — the two halves of a smooth gradient. There are no significant LH or HL units at all, which is itself informative: this surface is smooth, so no unit is unlike its surroundings.

On real data the spatial outliers (LH and HL) are often the most interesting rows in the table. A low-value unit surrounded by high-value ones is precisely the anomaly a global statistic averages away.

One caution: this runs a hypothesis test for every unit, so with 488 units and α = 0.05 you would expect about 24 false positives by chance — close to the 24 spatial outliers found above. The multiple-comparisons correction matters, and Moran_Local exposes p_z_sim and FDR-adjusted options for it.

Negative E[I] and why weights matter

Why E[I] is negative

Under the null hypothesis of no spatial pattern, the expected value of Moran's I is:

E[I] = −1 / (n − 1)

Not zero. The reason is that I is computed from deviations about the sample mean, and each unit's own value contributes to that mean — so the deviations are very slightly negatively correlated with each other by construction.

For 488 units, E[I] = −0.0021 and the distinction is cosmetic. For 20 units it is −0.053, and an observed I of −0.04 is positive autocorrelation relative to the null, despite being a negative number. Always compare against moran.EI, never against zero.

Why the weights are not a technicality

The definition of "neighbour" is the model of the spatial process. Contiguity says influence stops at a shared boundary; k-nearest-neighbours says every unit has the same number of influences regardless of geography; a distance band says influence has a range in metres.

The checkerboard result shows what happens when the definition and the process operate at different scales. And there is a second, quieter problem: contiguity weights on administrative units inherit all the arbitrariness of those units. Redraw the boundaries and the adjacency graph changes, and so does I.

The practical discipline: choose the definition from the process, test two or three alternatives, and report both the choice and the spread.

Why autocorrelation invalidates ordinary statistics

Independent observations each contribute one unit of information. Correlated observations do not — if you know a unit's neighbours, you already know most of what that unit will say.

The result is an effective sample size smaller than n. Standard errors computed as if all n observations were independent are too small, confidence intervals are too narrow, and p-values are too optimistic. It is not a small effect: with I near 0.9, most of the apparent information in 400 units is redundant.

This is why spatial regression exists. If a regression's residuals show significant Moran's I, the model is missing spatial structure, and the fix is a spatial lag or spatial error specification rather than tighter confidence intervals.

Why a global statistic can hide everything interesting

Moran's I is one number for the whole map. A study area with a strong hotspot in one corner and strong dispersion elsewhere can average to an I near zero — "no spatial pattern" — while containing two very strong ones.

That is the entire motivation for local indicators (LISA) and Getis-Ord Gi*: compute the statistic per unit, and map where the pattern actually is. Use the global statistic to answer "is there structure here at all", and local statistics to answer "where".

Islands, rates, residuals and Geary's C

• Islands — units with no neighbours — get zero weight and are effectively dropped, with a warning. See spatial weights warn about islands.
• use_index=True is required by current libpysal versions when building weights from a GeoDataFrame; without it you get a deprecation warning and index-alignment risk.
• Row-standardisation changes I slightly and changes the interpretation from "sum of neighbours" to "mean of neighbours". Row-standardise for irregular units.
• 999 permutations gives a minimum p of 0.001. For stricter thresholds use more permutations.
• Rates need care. Moran's I on a rate with wildly varying denominators picks up small-population noise. Use esda.moran.Moran_Rate, which applies an Empirical Bayes correction.
• Moran's I on regression residuals is the standard diagnostic for whether a model needs a spatial specification.
• Local Moran runs n tests. Correct for multiple comparisons, or expect roughly 5% false positives.
• Geary's C is an alternative that is more sensitive to local differences; the two often disagree at the margins, which is informative.

Internal links

• How to calculate Moran's I in Python — the practical implementation
• Spatial weights warn about islands or Moran's I returns nan — the commonest failure
• How to find hotspots with Getis-Ord Gi* in Python — locating what the global statistic averages away
• The modifiable areal unit problem explained — why the zones change the answer too
• How to bin points into hexagons in Python — a neighbour structure with no diagonal ambiguity
• Spatial join cardinality explained — the adjacency relationships underneath
• Kernel density explained — the continuous alternative to zone-based statistics
• How to aggregate spatial data by region in GeoPandas — producing the zone values in the first place
• A species distribution model scores well and predicts badly — autocorrelation inflating a model score

FAQ

What does Moran's I actually measure?

The correlation between each unit's value and the weighted average of its neighbours' values. Near +1 means similar values cluster; near −1 means they alternate.

Why is the expected value not zero?

Because I is computed from deviations about the sample mean. E[I] = −1/(n−1), which is negligible at n=488 and substantial at n=20. Always compare against moran.EI.

Which spatial weights should I use?

Choose from the process — contiguity for administrative units where boundaries matter, k-nearest for uneven unit sizes, a distance band when influence has a physical range. Test two or three and report the spread.

Why did my Moran's I change completely with different weights?

Because the pattern operates at the scale of the neighbourhood definition. A checkerboard gives I = −1.000 with rook weights and −0.039 with queen. If your result is that sensitive, the weights are the finding.

Should I use the analytical or permutation p-value?

The permutation one (p_sim). It makes no normality assumption. With 999 permutations the minimum reportable p is 0.001.

My Moran's I is significant. Now what?

If it was on raw values, map the local statistics to find where. If it was on regression residuals, your model is missing spatial structure and needs a spatial lag or error specification.

Can Moran's I be used on point data?

Not directly — it needs values attached to units with a neighbour structure. Aggregate the points to zones or a hex grid first, accepting the aggregation caveats.