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.

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.

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.

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