Getis-Ord Gi* Hot Spot Analysis in Python

TL;DR: Build a binary spatial weight matrix, then call esda.getisord.G_Local(y, w, star=True, permutations=999, seed=42). Read gi.Zs for the z-score (positive = hot spot, negative = cold spot) and gi.p_sim for the permutation p-value. Classify a unit as a hot or cold spot when p_sim clears your confidence threshold. Never row-standardise the weights for Gi*.

Why This Matters #

Getis-Ord $G_i^*$ turns a column of numbers on a map into a defensible answer to “where are the significant clusters?” It is the workhorse behind disease-cluster detection, crime hot spot policing, and pollution intensification maps. This walkthrough is the hands-on companion to the hot spot analysis guide, which develops the $G_i^*$ formula, permutation inference, and multiple-comparison correction in full. Both sit within the broader Core Concepts of Spatial Statistics & Geostatistics framework. Here we focus strictly on running the code end to end.

Environment #

pip install "geopandas>=1.0" "libpysal>=4.9" "esda>=2.5" "numpy>=1.23" "matplotlib>=3.7"

import geopandas as gpd
import numpy as np
import libpysal
from esda.getisord import G_Local
from shapely.geometry import box
import matplotlib.pyplot as plt

Step-by-Step Implementation #

Step 1 — Build or load a projected dataset #

# Synthetic 12x12 grid of polygon cells — replace with your own file:
# gdf = gpd.read_file("your_data.gpkg").to_crs("EPSG:26918")
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")
gdf = gdf.reset_index(drop=True)

The CRS must be projected (metres). Distance-band neighbourhoods computed on latitude/longitude are distorted, and reset_index(drop=True) guarantees the weight rows line up with the attribute rows.

Step 2 — Attach an attribute with an embedded hot spot #

rng = np.random.default_rng(42)
gdf["rate"] = rng.normal(loc=10.0, scale=2.0, size=len(gdf))

# Inject a high-value cluster in the lower-left quadrant so a hot spot exists
centroids = gdf.geometry.centroid
hot = (centroids.x < 4) & (centroids.y < 4)
gdf.loc[hot, "rate"] += 8.0

The attribute column must be numeric, free of NaN, and have non-zero variance. Constant columns make the z-score undefined.

Step 3 — Build BINARY spatial weights #

# Fixed distance band. binary=True is mandatory for Getis-Ord.
w = libpysal.weights.DistanceBand.from_dataframe(
    gdf, threshold=2.0, binary=True, silence_warnings=True
)
# Do NOT set w.transform = "R" — Gi* needs binary weights, not row-standardised.

DistanceBand returns binary weights by default. For areal polygons you can instead use libpysal.weights.Queen.from_dataframe(gdf), which is also binary. Row-standardising here would break the statistic.

Step 4 — Compute Gi* with star=True #

y = gdf["rate"].values.astype(float)
gi = G_Local(y, w, star=True, permutations=999, seed=42)

# gi.Zs    -> z-scores (one per unit)
# gi.p_sim -> permutation pseudo p-values
# gi.Gs    -> raw Gi* statistic before standardisation
gdf["z"] = gi.Zs
gdf["p"] = gi.p_sim

star=True adds the self-weight, giving $G_i^*$ rather than plain $G_i$. The explicit seed makes the permutation p-values reproducible.

Step 5 — Classify hot and cold spots at 90/95/99% #

def confidence_band(z, p):
    az = abs(z)
    if p <= 0.01 and az >= 2.58:
        band = 99
    elif p <= 0.05 and az >= 1.96:
        band = 95
    elif p <= 0.10 and az >= 1.65:
        band = 90
    else:
        return "Not Significant"
    return f"{'Hot' if z > 0 else 'Cold'} Spot {band}%"

gdf["gi_class"] = [confidence_band(z, p) for z, p in zip(gdf["z"], gdf["p"])]
print(gdf["gi_class"].value_counts())

Positive z-scores that clear a threshold become hot spots; negative z-scores become cold spots; everything else is not significant.

Step 6 — Map the result #

fig, ax = plt.subplots(figsize=(7, 7))
gdf.plot(column="z", cmap="RdBu_r", legend=True, ax=ax,
         edgecolor="0.8", linewidth=0.3, vmin=-4, vmax=4)
ax.set_title("Getis-Ord Gi* z-scores (red = hot, blue = cold)")
ax.set_axis_off()
fig.savefig("gi_star_hotspots.png", dpi=150, bbox_inches="tight")

A diverging colour map centred on zero is the correct choice: red for positive z (hot), blue for negative z (cold), white for units near zero. Fixing vmin and vmax symmetrically (here $\pm 4$) keeps the white midpoint anchored at $z = 0$ so the two directions stay visually comparable. For a publication map, colour by the gi_class categories instead of the raw z-score so that only statistically significant units are shaded and non-significant units read as neutral.

Interpreting the Output #

• gi.Zs — the signed z-score. Positive marks a hot spot, negative a cold spot; the injected lower-left cluster should show the largest positive values.
• gi.p_sim — the permutation pseudo p-value. With 999 permutations the smallest attainable value is 0.001. Only units below your threshold are genuine spots.
• gi.Gs — the raw, unstandardised statistic. Useful for diagnostics but not for classification; always classify on Zs and p_sim.
• Value counts — a healthy result shows spatially contiguous hot and cold blocks. A map of scattered isolated spots usually means the distance band is too tight or correction was skipped.

For the synthetic grid above, the injected lower-left cluster should surface as a compact patch of 99% hot spots, the far corner should trend cold or neutral, and the bulk of the map should read as not significant. If instead nearly every cell is flagged, the weights were row-standardised or no multiple-comparison correction was applied — both silently inflate the significant count. Compare the raw z magnitudes against the confidence bands: a cell at $z = 3.1$ belongs to a stronger, more coherent cluster than one at $z = 2.0$, even though both may clear the 95% threshold.

Critical Best Practices #

Keep the weights binary #

G_Local expects binary weights because $G_i^*$ sums neighbour values. Passing row-standardised weights silently changes the statistic’s meaning. Always build with binary=True and never set w.transform = "R". Weight-matrix construction is covered fully in spatial weight matrices.

Set a permutation seed #

Without a fixed seed, the permutation null — and therefore every p-value and every classified spot — changes on each run. Pass seed=42 (or any fixed integer) so results are reproducible in reports and tests.

Correct for multiple comparisons #

One test per unit means false positives accumulate. For maps beyond a few hundred units, apply a false discovery rate correction to p_sim before classifying, rather than a raw 0.05 cut. The parent hot spot analysis guide gives a Benjamini-Hochberg routine.

Choose the distance band deliberately #

The threshold encodes your hypothesis about interaction range. Too small and neighbourhoods collapse to single units (no hot spots); too large and everything blurs into one cluster. Inspect the nearest-neighbour distance distribution before fixing the band, and set the threshold above the maximum nearest-neighbour distance so no unit is left isolated:

from libpysal.weights import KNN
knn1 = KNN.from_dataframe(gdf, k=1)
nn_dist = np.array([knn1.neighbors[i][0] for i in range(len(gdf))])
# a safe lower bound for the band is the max 1-NN distance
print("min band to avoid islands:", round(float(nn_dist.max()), 2))

If the map still looks unstable, run Gi* at two or three bands and keep the hot spots that persist across scales — genuine clusters are robust to reasonable changes in the neighbourhood radius.

Troubleshooting #

Next Steps #

Return to the hot spot analysis guide for the mathematics and FDR correction, then compare against local correlation by working through the LISA workflow in how to calculate Moran’s I in PySAL.

← Back to Hot Spot Analysis with Getis-Ord Statistics

Related

• Hot Spot Analysis with Getis-Ord Statistics — the theory, permutation inference, and multiple-comparison correction behind this code
• Spatial Weight Matrices — building the binary weights Gi* requires
• How to Calculate Moran’s I in PySAL — the global-autocorrelation diagnostic to run first