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.34548422645435506 1.0107837784795915e-36

• 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

np.random.seed(0)                          # the permutation test draws at random: fix it to repeat
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.3455   E[I] = -0.0021
z = +13.48   p = 0.0010   n = 488

Five things every reported Moran's I needs:

Valid values, use_index, neighbours, row standardisation

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 0.0, max 174.37517954610743

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=488, mean neighbours 5.73, 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

Compare them on the same values. On the Leeds crime rates all three agree that the pattern is clustered, but not on how strongly:

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
}
for name, build in SCHEMES.items():
    wi = build(zones)
    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
  knn8      mean nb  8.0   I = +0.2786   p = 0.0010

Queen and rook have the same mean neighbour count here because real boundaries almost never meet at a single point; KNN, which ignores boundaries, gives a clearly lower I. On a pattern operating at the neighbourhood scale the schemes do not agree at all. A checkerboard is the cleanest example — a constructed 20 × 20 grid of alternating values, not data:

from shapely.geometry import box

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)])
for name, build, label in (("rook", Rook.from_dataframe, "4 edge neighbours"),
                           ("queen", Queen.from_dataframe, "8 neighbours")):
    wc = build(cells, use_index=True)
    wc.transform = "r"
    np.random.seed(0)
    mc = Moran(cells["value"].to_numpy(dtype="float64"), wc, permutations=999)
    print(f"  {name:<5} ({label:<17})  I = {mc.I:+.4f}   p = {mc.p_sim:.4f}")

  rook  (4 edge neighbours)  I = -1.0000   p = 0.0010
  queen (8 neighbours     )  I = -0.0393   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

np.random.seed(0)
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.

Reportable Moran's I with weights sensitivity

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 — 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]}")

    np.random.seed(0)
    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, "rate")

  queen n=488 mean_nb= 5.7 | I=+0.3455 E[I]=-0.0021 z= +13.48 p=0.0010 -> clustered

The verdict compares I against E[I], not against zero. At n=488 the difference is 0.0021 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, "rate")

  queen n=488 mean_nb= 5.7 | I=+0.3455 E[I]=-0.0021 z= +13.48 p=0.0010 -> clustered
  rook  n=488 mean_nb= 5.7 | I=+0.3570 E[I]=-0.0021 z= +13.77 p=0.0010 -> clustered
  knn   n=488 mean_nb= 4.0 | I=+0.2866 E[I]=-0.0021 z=  +9.96 p=0.0010 -> clustered
  knn   n=488 mean_nb= 8.0 | I=+0.2786 E[I]=-0.0021 z= +13.57 p=0.0010 -> clustered

weights  mean_nb      I     z     p
  queen      5.7 0.3455 13.48 0.001
   rook      5.7 0.3570 13.77 0.001
knn k=4      4.0 0.2866  9.96 0.001
knn k=8      8.0 0.2786 13.57 0.001

I spread across weights: 0.0784 (WEIGHTS-SENSITIVE — report all of them)

On Leeds the spread is 0.078 — KNN against the contiguity schemes — and the function flags it. On the checkerboard in step 3 the queen–rook gap alone is 0.96. The 0.05 threshold is arbitrary; the value of the line is that it forces the comparison to happen at all, and here it says to report more than one number.

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, "rate", w)
result.plot(column="cluster", categorical=True, legend=True,
            cmap="Set1", edgecolor="white", linewidth=0.2)

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

488 tests at α=0.05: expect ~24 false positives

59 low-low and 46 high-high units, and 24 spatial outliers: 15 low-rate LSOAs surrounded by high-rate ones (LH) and 9 the other way round (HL). Those two rows are usually the interesting ones — a low-crime area inside a high-crime cluster is exactly the anomaly the global I of 0.35 averages away.

The false-positive line matters. Running 488 tests at α = 0.05 produces about 24 significant results from noise alone — as many as the outliers found here — so treat isolated single-unit clusters with suspicion and contiguous blocks as real.

Misaligned weights and permutation p-values

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.

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
np.random.seed(0)
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.

Random seeds, Moran_Rate and LISA quadrants

• 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.

Internal links

• Spatial autocorrelation explained — what the statistic measures and why E[I] is not zero
• Spatial weights warn about islands or Moran's I returns nan — the failures this validation catches
• How to find hotspots with Getis-Ord Gi* in Python — locating what the global statistic averages away
• The modifiable areal unit problem explained — the zones the weights are built on
• How to bin points into hexagons in Python — turning points into units with neighbours
• How to aggregate spatial data by region in GeoPandas — producing the values
• How to handle missing and null values in spatial datasets — before the NaN reaches the statistic
• Choropleth classification explained — mapping the LISA categories

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.