βš™οΈ How-To Spatial statistics Joins & relationships Spatial Statistics Β· step 11 of 13 Updated

How to Calculate Moran's I in Python

Problem statement

You have values on a map and you need a number for "is this clustered?". Moran's I is that number, and the four-line version has three problems in it:

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

w = Queen.from_dataframe(zones)
moran = Moran(zones["rate"], w)
print(moran.I, moran.p_norm)
FutureWarning: `use_index` defaults to False but will default to True in future.
0.7876733636730491 0.0
  • A deprecation warning about index alignment that is not cosmetic β€” get it wrong and your values attach to the wrong polygons.
  • p_norm assumes normality, which spatial data rarely satisfies.
  • No w.transform, so a zone with twelve neighbours counts three times a zone with four.

And if the data contains a single NaN, moran.I comes back as nan while the p-value still looks significant.

Quick answer 

import numpy as np
from esda.moran import Moran
from libpysal.weights import Queen

values = zones["rate"].to_numpy(dtype="float64")
assert np.isfinite(values).all(), "non-finite values make Moran's I nan"
assert np.var(values) > 0, "zero variance makes Moran's I undefined"

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

moran = Moran(values, w, permutations=999)
print(f"I = {moran.I:+.4f}   E[I] = {moran.EI:+.4f}")
print(f"z = {moran.z_sim:+.2f}   p = {moran.p_sim:.4f}   n = {w.n - len(w.islands)}")
I = +0.9046   E[I] = -0.0025
z = +34.52   p = 0.0010   n = 400

Five things every reported Moran's I needs:

Report Why
I and E[I] E[I] is βˆ’1/(nβˆ’1), not zero
the weights scheme queen, rook and KNN give different answers
p_sim, not p_norm permutation inference makes no distributional assumption
the permutation count 999 permutations means the minimum p is 0.001
n actually used islands are silently excluded
Five steps to a reportable Moran's I: validate the values, build weights, row-standardise, run permutations, report with the weights scheme.
Two of the five steps are validation. They are the two most often skipped.

Step-by-step solution

1. Validate the values before anything else 

values = zones["rate"].to_numpy(dtype="float64")
bad = ~np.isfinite(values)
if bad.any():
    raise ValueError(f"{bad.sum()} non-finite values at {zones.index[bad].tolist()[:5]}")

Moran's I sums over every unit, so one NaN or inf makes the whole statistic nan. What makes this dangerous is that p_sim still returns 0.001 β€” see spatial weights warn about islands or Moran's I returns nan for why.

Also check for suppression sentinels, which are finite and therefore pass the test above:

print(f"min {values.min()}, max {values.max()}")
min -99999.0, max 84.2

2. Build the weights with use_index=True 

w = Queen.from_dataframe(zones, use_index=True)
print(f"n={w.n}, mean neighbours {w.mean_neighbors:.2f}, islands {len(w.islands)}")
n=400, mean neighbours 7.41, islands 0

use_index=True keys the weights on the GeoDataFrame's index rather than positional order. Without it, any reindexing between building the weights and computing the statistic silently misaligns values and polygons β€” and the result looks entirely plausible.

3. Choose the neighbour definition deliberately 

from libpysal.weights import KNN, Queen, Rook

SCHEMES = {
    "queen": lambda g: Queen.from_dataframe(g, use_index=True),   # shares a boundary or corner
    "rook": lambda g: Rook.from_dataframe(g, use_index=True),     # shares an edge only
    "knn8": lambda g: KNN.from_dataframe(g, k=8),                 # 8 nearest, regardless of touching
}

On a smooth surface these agree closely:

  queen     mean nb  7.4   I = +0.9046   p = 0.0010
  rook      mean nb  3.8   I = +0.9152   p = 0.0010
  knn8      mean nb  8.0   I = +0.9005   p = 0.0010

On a pattern operating at the neighbourhood scale they do not agree at all:

  A checkerboard
  rook  (4 edge neighbours)   I = -0.9980   p = 0.0010
  queen (8 neighbours)        I = -0.0386   p = 0.0700

Run two or three and report the spread. A wide spread means the weights are the finding.

4. Row-standardise 

w.transform = "r"

Each unit's weights now sum to 1, so the spatial lag is the mean of its neighbours rather than the sum. Without it, units with many neighbours dominate purely because of how the boundaries fall β€” which for irregular administrative units is a large and arbitrary effect.

5. Use permutation inference 

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

p_sim shuffles the values across the units 999 times and counts how often a random arrangement is as extreme as yours. p = 0.0010 means "never in 999 tries" β€” the smallest value 999 permutations can produce, not a precise probability.

Moran's I on a smooth surface agreeing across queen, rook and KNN weights, and disagreeing completely on a checkerboard.
When the schemes agree, the result is robust. When they do not, the neighbour definition is the result.

Code examples

Example 1 β€” a function that returns something reportable 

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

SCHEMES = {
    "queen": lambda g, k: Queen.from_dataframe(g, use_index=True),
    "rook": lambda g, k: Rook.from_dataframe(g, use_index=True),
    "knn": lambda g, k: KNN.from_dataframe(g, k=k),
}


def morans_i(zones, column, *, scheme="queen", k=8, permutations=999, transform="r"):
    """Global Moran's I with every choice made explicit and validated."""
    values = zones[column].to_numpy(dtype="float64")

    bad = ~np.isfinite(values)
    if bad.any():
        raise ValueError(f"{bad.sum()} non-finite value(s) in {column!r} at index "
                         f"{zones.index[bad].tolist()[:5]}")
    if np.var(values) == 0:
        raise ValueError(f"{column!r} has zero variance β€” Moran's I is undefined")
    if values.min() < 0 and column.endswith(("rate", "count", "density")):
        print(f"  note: {column!r} has negative values (min {values.min():,.0f}) β€” "
              f"check for suppression sentinels")

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        w = SCHEMES[scheme](zones, k)
    w.transform = transform

    if w.islands:
        print(f"  {len(w.islands)} island(s) excluded: {list(w.islands)[:5]}")

    moran = Moran(values, w, permutations=permutations)
    if not np.isfinite(moran.I):
        raise ValueError(f"Moran's I is {moran.I} despite finite inputs")

    used = w.n - len(w.islands)
    verdict = ("clustered" if moran.I > moran.EI else "dispersed") \
        if moran.p_sim <= 0.05 else "no significant pattern"

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


moran, w = morans_i(zones, "value")
  queen n=400 mean_nb= 7.4 | I=+0.9046 E[I]=-0.0025 z= +34.52 p=0.0010 -> clustered

The verdict compares I against E[I], not against zero. At n=400 the difference is 0.0025 and cosmetic; at n=20 it is 0.053, and an I of βˆ’0.04 is genuinely positive autocorrelation relative to the null.

Example 2 β€” reporting the weights sensitivity 

def sensitivity(zones, column, schemes=(("queen", 8), ("rook", 8), ("knn", 4), ("knn", 8))):
    rows = []
    for scheme, k in schemes:
        moran, w = morans_i(zones, column, scheme=scheme, k=k)
        rows.append({
            "weights": f"{scheme}" + (f" k={k}" if scheme == "knn" else ""),
            "mean_nb": round(w.mean_neighbors, 1),
            "I": round(moran.I, 4),
            "z": round(moran.z_sim, 2),
            "p": moran.p_sim,
        })
    frame = pd.DataFrame(rows)
    spread = frame["I"].max() - frame["I"].min()
    print(f"\n{frame.to_string(index=False)}")
    print(f"\nI spread across weights: {spread:.4f} "
          f"({'robust' if spread < 0.05 else 'WEIGHTS-SENSITIVE β€” report all of them'})")
    return frame


sensitivity(zones, "value")
  weights  mean_nb       I      z      p
    queen      7.4  0.9046  34.52  0.001
     rook      3.8  0.9152  34.19  0.001
  knn k=4      4.0  0.9167  36.02  0.001
  knn k=8      8.0  0.9005  40.13  0.001

I spread across weights: 0.0162 (robust)

Run the same function on a checkerboard and the spread is 0.96, and the label changes to WEIGHTS-SENSITIVE. That threshold is arbitrary β€” the value of the line is that it forces the comparison to happen at all.

Example 3 β€” local Moran, to find out where 

from esda.moran import Moran_Local


def local_moran(zones, column, w, *, permutations=999, alpha=0.05):
    """LISA: per-unit cluster type, with the non-significant units labelled as such."""
    values = zones[column].to_numpy(dtype="float64")
    np.random.seed(42)                        # permutation results move between runs
    lisa = Moran_Local(values, 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]

    counts = out["cluster"].value_counts()
    expected_false = int(alpha * len(zones))
    print(counts.to_string())
    print(f"\n{len(zones)} tests at Ξ±={alpha}: expect ~{expected_false} false positives")
    return out


result = local_moran(zones, "value", w)
result.plot(column="cluster", categorical=True, legend=True,
            cmap="Set1", edgecolor="white", linewidth=0.2)
cluster
not significant    173
LL                 116
HH                 111

400 tests at Ξ±=0.05: expect ~20 false positives

Two large concentrations and no spatial outliers at all β€” no LH or HL units. That is the correct answer for a smooth gradient: nothing is unlike its surroundings.

On real, messier data the LH and HL rows are usually the interesting ones. A low-value unit surrounded by high-value neighbours is exactly the anomaly the global statistic averages into a single number.

The false-positive line matters. Running 400 tests at Ξ± = 0.05 produces about 20 significant results from noise alone, so treat isolated single-unit clusters with suspicion and contiguous blocks as real.

Explanation

Why use_index=True is not a cosmetic warning

Queen.from_dataframe builds a neighbour graph keyed either on positional order or on the DataFrame index. If it uses positions, and you later filter, sort or reindex the frame before computing the statistic, values and polygons come apart.

The failure is silent. Moran's I on shuffled values is a valid number β€” it will be near E[I], which reads as "no spatial pattern", so a genuinely clustered dataset can report as random. Setting use_index=True keys everything on the index and removes the class of error.

Why row-standardisation matters most for irregular units

On a regular grid every interior cell has the same number of neighbours, so standardising changes almost nothing.

On administrative units, neighbour counts vary from 1 to 15 or more. Without standardisation the spatial lag is a sum, so a unit with 15 neighbours contributes fifteen times as much to the statistic as a unit with one β€” reflecting nothing but the shape of the boundaries. Row-standardising converts the lag to a mean and removes the artefact.

The trade-off: row-standardised weights are no longer symmetric, which matters for some spatial regression specifications but not for Moran's I.

Why permutation p-values are the default

The analytical p-value comes from a normal approximation to the sampling distribution of I. That approximation is good when n is large and the values are roughly normal, and unreliable when either fails β€” which for counts, rates and skewed spatial variables is most of the time.

Conditional permutation makes no assumption: shuffle the values across the units, recompute I, repeat. The observed I is compared against that empirical distribution. It costs 999 recomputations, which is milliseconds, and it is why p_sim exists alongside p_norm.

Note the resolution limit. With 999 permutations the smallest possible p is 1/1000 = 0.001. If you need to claim p < 0.0001, use 9,999 permutations.

Local Moran assigning each unit to HH, LL, LH, HL or not significant, with a smooth surface producing only HH and LL.
A smooth gradient produces only HH and LL. Spatial outliers β€” LH and HL β€” are what a global statistic hides.

Why to compute Moran's I on regression residuals

The most valuable use of this statistic is not on raw values at all. Fit a model, take the residuals, and test them:

residuals = model.resid
moran_resid = Moran(residuals, w, permutations=999)
print(f"residual I = {moran_resid.I:+.4f}, p = {moran_resid.p_sim:.4f}")

A significant result means neighbouring observations have similar errors β€” the model is missing spatial structure. Its standard errors are then too small and its p-values too optimistic, because the residuals are not independent.

The fix is a spatial specification (spatial lag or spatial error), not a tighter confidence interval. This diagnostic is the standard reason spatial econometrics exists.

Edge cases or notes

  • use_index=True is required on current libpysal versions and prevents silent misalignment.
  • Set np.random.seed() before permutation tests if the exact numbers must be reproducible; results move by a few units otherwise.
  • Moran_Rate applies an Empirical Bayes correction and is the right tool for rates with varying denominators β€” a rate over 12 people is mostly noise.
  • Islands are excluded silently. Report w.n - len(w.islands).
  • Local Moran runs n tests. Expect roughly Ξ±Β·n false positives, and prefer contiguous blocks over isolated units.
  • lisa.q is 1–4 for HH, LH, LL, HL and is populated for every unit regardless of significance β€” mask it with p_sim before using it.
  • Geary's C is more sensitive to local differences than Moran's I. Where they disagree, look at the local statistics.
  • Moran's I needs areal units with neighbours. For raw points, aggregate to zones or a hex grid first.

FAQ

What library should I use?

libpysal for the weights and esda for the statistic. Both are part of the PySAL family and are the standard tools.

Why is my Moran's I nan?

A non-finite value in the data, or zero variance. Islands do not cause it. Check np.isfinite(values).all() before running, because p_sim will still return a significant-looking number.

Should I use p_norm or p_sim?

p_sim. The analytical p-value assumes normality, which spatial counts and rates rarely satisfy. Permutation inference makes no such assumption and costs milliseconds.

What does w.transform = "r" do?

Row-standardises the weights so each unit's sum to 1, making the spatial lag a mean rather than a sum. Important for irregular units where neighbour counts vary widely.

Which weights scheme should I pick?

Whichever matches the process. Run two or three and report the spread β€” if the spread is small the result is robust, and if it is large the weights are the finding.

How do I find where the clusters are?

Moran_Local (LISA), which classifies each unit as HH, LL, LH, HL or not significant. Correct for multiple comparisons or expect about 5% false positives.

What does a significant Moran's I on residuals mean?

Your model is missing spatial structure, so its standard errors are too small. Move to a spatial lag or spatial error specification.