Geostatistical Modeling Spatial Statistics with Python

Computing Local Moran's I (LISA) in Python

TL;DR: Call esda.moran.Moran_Local(y, w, permutations=999) on a row-standardized weights matrix. The result exposes Is (per-unit local statistics), q (quadrant codes 1–4 = HH/LH/LL/HL), and p_sim (per-unit pseudo p-values). Keep only units where p_sim < 0.05, map their quadrant to a spatial-cluster label, and plot the result as a LISA cluster map.

Why This Matters

A significant global Moran’s I tells you spatial autocorrelation exists somewhere in your study area, but not where. Local Indicators of Spatial Association (LISA) decompose that single global number into one statistic per location, letting you pinpoint hot spots, cold spots, and spatial outliers on a map. This is the step every applied analyst needs after the global test: it turns “there is clustering” into “these specific tracts form a high-high cluster.” Local Moran’s I is a core member of the spatial autocorrelation metrics family within the broader Core Concepts of Spatial Statistics & Geostatistics guide, and it is the natural follow-up to a confirmed global result.

The local statistic for location ii is:

Ii=(xi−xˉ)m2∑jwij(xj−xˉ),m2=∑i(xi−xˉ)2nI_i = \frac{(x_i - \bar{x})}{m_2} \sum_j w_{ij}(x_j - \bar{x}), \qquad m_2 = \frac{\sum_i (x_i - \bar{x})^2}{n}

The sum of all IiI_i is proportional to the global Moran’s I, which is why LISA is a genuine decomposition rather than a separate statistic. The sign of IiI_i combined with whether xix_i is above or below the mean gives the four-way quadrant classification that drives the LISA cluster map.

Where esda's quadrant codes come from Left: a Moran scatterplot with each of the 144 grid cells plotted as its own standardized value against the mean of its neighbours, split into four quadrants labelled q = 1 High-High, q = 2 Low-High, q = 3 Low-Low and q = 4 High-Low; a dashed regression line through the origin has slope 0.58, the global Moran's I. Right: four cards restating each quadrant code, its cluster or outlier meaning, and the sign of the local statistic. The quadrant codes are a position on the Moran scatterplot one point per grid cell – its own standardized value against its neighbours’ mean; the dashed slope is the global Moran’s I, 0.58 q = 2 Low-High q = 1 High-High q = 3 Low-Low q = 4 High-Low z(i): the unit’s own value, standardized (0 = study mean) Wz(i): neighbour mean q = 1 → High-High – cluster core above the mean, neighbours above → I(i) > 0 q = 2 → Low-High – spatial outlier below the mean, neighbours above → I(i) < 0 q = 3 → Low-Low – cluster core below the mean, neighbours below → I(i) > 0 q = 4 → High-Low – spatial outlier above the mean, neighbours below → I(i) < 0 faded points have p_sim > 0.05 – they still carry a q code, which is why masking matters

Environment and Version Pinning

bash
pip install "geopandas>=1.0" "libpysal>=4.9.0" "esda>=2.5.0" "numpy>=1.23" "matplotlib>=3.7"
python
import geopandas as gpd
import numpy as np
import libpysal
from esda.moran import Moran_Local
import matplotlib.pyplot as plt
from shapely.geometry import box

Step-by-Step Implementation

Step 1 — Load data and attach the attribute

python
# Synthetic 12x12 grid with a diagonal gradient so real clusters appear —
# replace with gdf = gpd.read_file("your_data.gpkg")
cells = [box(x, y, x + 1, y + 1) for y in range(12) for x in range(12)]
gdf = gpd.GeoDataFrame(geometry=cells, crs="EPSG:32618").reset_index(drop=True)

rng = np.random.default_rng(11)
cx = gdf.geometry.centroid.x.values
cy = gdf.geometry.centroid.y.values
gdf["pm25"] = 6.0 + 0.5 * (cx + cy) + rng.normal(0, 2.0, size=len(gdf))

The attribute must be numeric with no NaN, and the geometry must sit in a projected CRS so contiguity is computed on true adjacency rather than distorted degrees.

Step 2 — Build and row-standardize the weights matrix

python
w = libpysal.weights.Queen.from_dataframe(gdf)
w.transform = "R"
y = gdf["pm25"].values

Row-standardization makes each unit’s spatial lag a weighted average of its neighbours, which is required for the quadrant logic to be comparable across units with different neighbour counts. For the trade-off between Queen and Rook adjacency here, see Rook vs Queen Contiguity Weights.

Step 3 — Run Moran_Local

python
lisa = Moran_Local(y, w, permutations=999, seed=42)

Moran_Local computes one local statistic per observation and runs a conditional permutation test for each: it holds unit ii fixed and reshuffles the remaining values among its neighbours to build a per-unit null distribution. seed makes the pseudo p-values reproducible.

What the conditional permutation test holds fixed Left: two three-by-three neighbourhoods. In the observed one the centre cell reads 16.2 and its eight neighbours average 15.27, giving a local Moran's I of 1.47; in a permuted one the same centre value 16.2 is kept while eight values drawn from elsewhere in the grid average 11.96, giving a local I of minus 0.03. Right: the null distribution built from 999 such draws, centred near zero, with the observed 1.47 marked in the far right tail where only two draws reach, so p_sim is 0.003. Moran_Local tests each unit against its own reshuffled neighbourhood 1. redraw the neighbourhood 999 times 2. read the observed I(i) against that null unit i (red) and 8 neighbours same i, neighbours redrawn 16.6 16.1 17.0 13.9 16.2 14.5 17.2 13.6 13.1 12.4 13.8 8.4 14.5 16.2 11.3 9.9 18.2 7.2 × 999 observed nbr mean 15.27, I(i) = 1.47 draw 1 of 999 nbr mean 11.96, I(i) = −0.03 the centre value never moves – only the eight neighbour values are redrawn from elsewhere in the grid, so the null asks whether this unit’s surroundings are unusual, not the unit alone study-area mean pm25 = 12.05, m2 = 9.07 −2 −1 0 1 2 observed I(i) = 1.47 only 2 of the 999 draws reach it local I(i) across 999 conditional permutations p_sim = (2 + 1) / (999 + 1) = 0.003 the smallest p_sim 999 permutations can resolve is 0.001

Step 4 — Extract the Is, quadrant codes, and p-values

python
gdf["local_I"] = lisa.Is        # local Moran's I per unit
gdf["quadrant"] = lisa.q        # 1=HH, 2=LH, 3=LL, 4=HL
gdf["p_sim"] = lisa.p_sim       # permutation pseudo p-value per unit

print(gdf[["local_I", "quadrant", "p_sim"]].head())

lisa.q is a NumPy integer array. Its coding is fixed by esda: 1 = High-High, 2 = Low-High, 3 = Low-Low, 4 = High-Low.

Step 5 — Classify quadrants and filter by significance

python
labels = {1: "High-High", 2: "Low-High", 3: "Low-Low", 4: "High-Low"}
alpha = 0.05

gdf["cluster"] = np.where(
    gdf["p_sim"] <= alpha,
    gdf["quadrant"].map(labels),
    "Not significant",
)

print(gdf["cluster"].value_counts())

Only units passing the p_sim <= alpha test keep their quadrant label; everything else collapses to "Not significant". This masking is what separates a real LISA cluster map from an over-interpreted one.

Step 6 — Plot the LISA cluster map

python
cluster_colors = {
    "High-High": "#c1272d",
    "Low-Low": "#0000a7",
    "Low-High": "#5ea3d1",
    "High-Low": "#f4a582",
    "Not significant": "#d9d9d9",
}
fig, ax = plt.subplots(figsize=(7, 7))
for label, color in cluster_colors.items():
    subset = gdf[gdf["cluster"] == label]
    if not subset.empty:
        subset.plot(ax=ax, color=color, edgecolor="white", linewidth=0.3, label=label)
ax.legend(loc="upper left", fontsize=8)
ax.set_axis_off()
ax.set_title("Local Moran's I (LISA) cluster map")
fig.savefig("lisa_clusters.png", dpi=150, bbox_inches="tight")

Interpreting the Output

The four significant categories carry distinct meaning:

Quadrant q code Meaning Interpretation
High-High 1 High value, high neighbours Hot spot — core of a spatial cluster
Low-Low 3 Low value, low neighbours Cold spot — core of a spatial cluster
Low-High 2 Low value, high neighbours Spatial outlier (negative local association)
High-Low 4 High value, low neighbours Spatial outlier (negative local association)

High-High and Low-Low are spatial clusters: the unit and its surroundings reinforce each other. Low-High and High-Low are spatial outliers: a unit that disagrees with its neighbourhood, which often flags data errors, boundary effects, or genuinely anomalous locations worth investigating. The magnitude of local_I ranks how strongly each unit contributes, but the cluster label — not the raw local_I — is what you communicate on the map.

Expect most units to be “Not significant.” That is normal and correct: LISA concentrates the eye on the minority of locations where local structure is statistically credible.

Critical Best Practices

Always mask by p_sim before mapping

lisa.q assigns a quadrant to every unit, significant or not. Plotting raw quadrants produces a map that looks fully clustered even under pure noise. Filter with p_sim <= alpha first; the un-masked quadrant array is only an intermediate.

Raw lisa.q against the significance-masked map Left: the 12 by 12 grid with every one of its 144 cells filled by its raw quadrant code, 56 High-High, 57 Low-Low, 16 Low-High and 15 High-Low, so the whole map appears clustered. Right: the same grid after masking on p_sim less than or equal to 0.05, where only 34 High-High, 31 Low-Low, 2 Low-High and 2 High-Low cells keep their colour and the remaining 75 are grey. A legend on the far right names the five fill colours. Every cell gets a quadrant; only 69 of the 144 earn a colour raw lisa.q – nothing masked masked where p_sim ≤ 0.05 56 HH, 57 LL, 16 LH, 15 HL – not one cell left grey pure noise would fill this map just as completely 34 HH, 31 LL, 2 LH, 2 HL – 75 cells not significant colour now survives only where the permutation test backs it cluster label High-High (q = 1) Low-Low (q = 3) Low-High (q = 2) High-Low (q = 4) Not significant same y, same weights, same lisa.q – only the p_sim mask differs

Correct for multiple comparisons

You run one test per location, so with hundreds of units the family-wise false-positive rate balloons. Apply a False Discovery Rate adjustment or tighten the per-test threshold to 0.01. A simple conservative option is to compare p_sim against alpha / n_tests.

Set a seed and enough permutations

Conditional permutation is stochastic. Pass seed= for reproducibility and raise permutations to 9999 for publication-grade stability, since the smallest resolvable pseudo p-value is 1/(permutations+1)1/(\text{permutations}+1).

Interpret outliers, do not discard them

Low-High and High-Low units are the most information-rich part of a LISA map. Before treating them as noise, check whether they represent measurement errors, edge units with truncated neighbourhoods, or real local anomalies — each demands a different response.

Keep global and local analysis consistent

Use the same attribute and the same weights object for both the global Moran’s I and the local decomposition. If the global result was not significant, be cautious reading local clusters, and cross-check the choice of statistic against Moran’s I vs Geary’s C: Which to Use.

Troubleshooting

Symptom Likely cause Fix
Entire map looks clustered Plotting lisa.q without significance masking Filter on p_sim <= alpha before assigning labels
Almost no significant units Too few permutations or very weak signal Raise permutations; verify the global Moran’s I is significant first
Quadrant labels look swapped Assumed custom coding Use the fixed esda coding: 1=HH, 2=LH, 3=LL, 4=HL
UserWarning: islands Disconnected units break the local lag Switch to KNN.from_dataframe(gdf, k=8) or remove islands
Results change every run No random seed set Pass seed= to Moran_Local
Many spurious clusters in large data No multiple-testing correction Apply FDR or use alpha=0.01

Next Steps

Confirm the global signal first with How to Calculate Moran’s I in PySAL, then decide whether a global product-moment or a difference-based index best fits your question via Moran’s I vs Geary’s C: Which to Use. Both return to the spatial autocorrelation metrics topic.

Frequently Asked Questions

What do the four LISA quadrants mean?

The quadrant code in lisa.q classifies each significant location by its own value and its neighbours’ average. High-High (1) is a hot spot surrounded by high values, Low-High (2) is a low outlier among high neighbours, Low-Low (3) is a cold spot surrounded by low values, and High-Low (4) is a high outlier among low neighbours. HH and LL are spatial clusters; HL and LH are spatial outliers.

Why must I filter LISA results by p_sim?

Every location receives a quadrant label whether or not its local statistic is significant, so an unfiltered map shows structure that may be noise. Mask locations where p_sim exceeds your threshold (commonly 0.05) and label only the significant units, otherwise you over-interpret random fluctuation as clustering.

Should I correct LISA p-values for multiple testing?

Yes, because you run one hypothesis test per location, the family-wise error rate inflates with nn. Apply a False Discovery Rate or Bonferroni-style adjustment, or use a stricter per-test threshold such as 0.01, to avoid flagging spurious clusters when the dataset is large.

← Back to Spatial Autocorrelation Metrics

Related

Related

How to Calculate Moran's I in PySAL Moran's I vs Geary's C: Which to Use