Geostatistical Modeling Spatial Statistics with Python

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 Gi∗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 Gi∗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

bash
pip install "geopandas>=1.0" "libpysal>=4.9" "esda>=2.5" "numpy>=1.23" "matplotlib>=3.7"
python
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

python
# 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

python
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

python
# 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.

Why Gi* needs binary weights On the left, a five-by-five patch of unit cells with a dashed circle of radius 2.0 drawn around the focal centroid: twelve neighbouring centroids fall inside and receive weight one, twelve fall outside and receive weight zero, and the focal cell adds its own self-weight because star=True. On the right, two stacked panels show the same thirteen values of about 18.0 summed with binary weights, giving 234.0 with a weight total of 13, and then with row-standardised weights of one thirteenth, giving 18.0 with a weight total of 1, which collapses the sum into a local mean. Gi* sums neighbour values — so the weights must stay binary Who is a neighbour at threshold = 2.0? i r = 2.0 w = 1 — centroid inside the band (12 cells) self-weight 1 — star=True adds cell i w = 0 — centroid beyond 2.0, excluded binary=True — every included weight is 1 Σ w·y = 13 × 18.0 = 234.0 Σ w = 13 13 cells of the injected hot quadrant, rate ≈ 18.0 Gi* asks: is that SUM bigger than 13 random draws? w.transform = 'R' — every weight becomes 1/13 Σ w·y = 13 × (18.0 ÷ 13) = 18.0 Σ w = 1 the sum has collapsed into a local mean every unit now reports Σ w = 1, and Gi* loses its scale Row-standardise for Moran's I; never for Getis-Ord Gi*.

Step 4 — Compute Gi* with star=True

python
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 Gi∗G_i^* rather than plain GiG_i. The explicit seed makes the permutation p-values reproducible.

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

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

How confidence_band classifies a unit Four stacked rows. The first three rows each hold one test: p_sim at most 0.01 with absolute z at least 2.58, then 0.05 with 1.96, then 0.10 with 1.65. Passing a test sends the unit to a hot spot label at that confidence when z is positive and a cold spot label when z is negative; failing drops it to the next row, and a unit that fails all three is labelled Not Significant and left unshaded. confidence_band(z, p): the first test that passes wins z > 0 → hot z < 0 → cold TEST 1 p_sim ≤ 0.01 and |z| ≥ 2.58 pass Hot Spot 99% e.g. z = +3.10, p = 0.001 Cold Spot 99% e.g. z = −2.95, p = 0.001 fail TEST 2 p_sim ≤ 0.05 and |z| ≥ 1.96 pass Hot Spot 95% e.g. z = +2.00, p = 0.012 Cold Spot 95% e.g. z = −2.10, p = 0.028 fail TEST 3 p_sim ≤ 0.10 and |z| ≥ 1.65 pass Hot Spot 90% e.g. z = +1.72, p = 0.078 Cold Spot 90% e.g. z = −1.80, p = 0.066 fail NO TEST PASSES p_sim > 0.10 or |z| < 1.65 Not Significant left unshaded — the sign of z carries no evidence here

Step 6 — Map the result

python
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 ±4\pm 4) keeps the white midpoint anchored at z=0z = 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.1z = 3.1 belongs to a stronger, more coherent cluster than one at z=2.0z = 2.0, even though both may clear the 95% threshold.

Three Gi* maps of the same 144 cells Three twelve-by-twelve grids shaded by Gi* class. The first, built with a 2.0 binary distance band, flags 31 of 144 cells: a compact nine-cell core of 99 percent hot spots in the lower-left quadrant fading through 95 and 90 percent, with a few cold cells in the far upper-right corner and the bulk of the map not significant. The second, with the band tightened to 1.0, flags only 11 scattered isolated cells with no coherent cluster. The third, with row-standardised weights and no false-discovery correction, floods 138 of 144 cells with hot or cold shading. Three Gi* maps of the same 144 cells — only one is trustworthy Healthy — band 2.0, binary 31 of 144 flagged; one compact core Band too tight — threshold 1.0 11 flagged, all isolated singletons Row-standardised, no FDR 138 of 144 flagged — not credible hot 99% hot 95% hot 90% not sig. cold 90% cold 95% cold 99%

Critical Best Practices

Keep the weights binary

G_Local expects binary weights because Gi∗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:

python
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

Symptom Likely cause Fix
Zs full of nan Constant attribute or island units with no neighbours Check y.std(); use KNN.from_dataframe(gdf, k=8) weights
Nearly every unit significant Row-standardised weights or no correction Rebuild with binary=True; apply an FDR threshold
No hot spots at all Distance band too small Increase threshold; inspect nearest-neighbour distances
Results change each run Permutation seed not set Pass seed= to G_Local
ValueError: shape mismatch Index not aligned with weights gdf.reset_index(drop=True) before building w
AttributeError: 'W' object has no attribute Legacy import path Use from esda.getisord import G_Local

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

Related

Emerging Hot Spot Analysis Over Time in Python Getis-Ord Gi* vs Local Moran's I for Hot Spots