βš™οΈ How-To Spatial statistics Maps & visualisation Spatial Statistics Β· step 13 of 13 Updated

How to Find Hotspots with Getis-Ord Gi* in Python

Problem statement

A heatmap shows where events are. It does not tell you whether any of those concentrations is more than chance would produce β€” kernel density finds peaks in uniformly random points too.

Getis-Ord Gi* is the statistic that answers the harder question: for each unit, is the local sum of values significantly higher (or lower) than you would expect if the values were shuffled at random across the map?

The standard recipe you will find in most tutorials looks like this:

from esda.getisord import G_Local

gi = G_Local(zones["value"].to_numpy(), w, star=True, permutations=999)
zones["hot"] = gi.Zs > 1.96          # "95% confidence"
print(zones["hot"].sum(), "hotspots")
0 hotspots

Zero β€” on data with Moran's I of 0.90 and a p-value of 0.001. The pattern is unambiguously clustered and this test found nothing at all.

The problem is Zs. On this data its analytical standard deviation is 3.1 times larger than the permutation-based one, so the z-scores are compressed into Β±1.83 and nothing ever clears 1.96.

Quick answer

Use the permutation p-value, and halve the threshold because it is one-sided:

import numpy as np
from esda.getisord import G_Local

gi = G_Local(zones["value"].to_numpy(), w, star=True, permutations=999)

z_sim = (gi.Gs - gi.EG_sim) / gi.seG_sim          # permutation-based z
significant = gi.p_sim <= 0.025                    # one-sided p, so Ξ±/2

zones["gi_z"] = z_sim
zones["hotspot"] = np.where(~significant, "not significant",
                            np.where(z_sim > 0, "hot", "cold"))
print(zones["hotspot"].value_counts().to_string())
hotspot
hot                101
cold               104
not significant    195
Field Use it? Why
Zs, p_norm no analytical variance is unreliable β€” found 0 hotspots on strongly clustered data
p_sim yes conditional permutation; one-sided, so compare against Ξ±/2
(Gs - EG_sim) / seG_sim yes permutation z-score, for the hot/cold direction and for mapping
Analytical Getis-Ord z-scores spanning only plus or minus 1.83 and finding zero hotspots, against permutation z-scores spanning plus or minus 5.2 and finding 205.
Same data, same weights. The analytical variance is over three times too large, so nothing ever reaches significance.

Step-by-step solution

1. Get your data into units with a neighbour structure

Gi* needs a value per areal unit and a definition of which units are neighbours. Points must be aggregated first β€” to administrative zones, a hex grid, or a square grid:

hexes = h3_bin(incidents, resolution=9)          # counts per cell

Everything from the modifiable areal unit problem applies: the resolution you choose changes what counts as a hotspot. Pick it from the process and say what you picked.

2. Build the weights, and think about what "nearby" means 

from libpysal.weights import Queen, DistanceBand, KNN

w = Queen.from_dataframe(hexes, use_index=True)
w.transform = "r"
print(f"n={w.n}, mean neighbours {w.mean_neighbors:.1f}, islands {len(w.islands)}")
n=400, mean neighbours 7.4, islands 0

For Gi* specifically, a distance band often models the question better than contiguity: a hotspot is usually "elevated over some spatial scale", and the band makes that scale explicit in metres rather than implicit in the shape of the zones.

w = DistanceBand.from_dataframe(hexes, threshold=800, binary=True)

Islands are fatal for a distance band β€” a unit with no neighbours inside the threshold has no local sum. Check w.islands and see spatial weights warn about islands.

3. Use star=True 

gi = G_Local(values, w, star=True, permutations=999)

star=True includes the unit itself in its own neighbourhood β€” that is what makes it Gi* rather than Gi. It is almost always what you want: a hotspot ought to include the hot unit, not just its surroundings.

With star=False a high-value unit surrounded by low ones scores low, which is a legitimate statistic answering a different question and rarely the one asked.

4. Ignore the analytical z-scores 

print(f"analytical Zs range  {gi.Zs.min():+.2f} .. {gi.Zs.max():+.2f}")
print(f"permutation z range  {z_sim.min():+.2f} .. {z_sim.max():+.2f}")
print(f"analytic sd / sim sd {(np.sqrt(gi.VGs) / gi.seG_sim).mean():.2f}")
analytical Zs range  -1.67 .. +1.83
permutation z range  -4.81 .. +5.17
analytic sd / sim sd 3.11

The analytical variance formula assumes a normal distribution of values. Real spatial data is skewed and autocorrelated, and here the formula overestimates the standard deviation by a factor of three, crushing every z-score toward zero.

The conditional permutation approach makes no such assumption: it shuffles the observed values across the units and measures the actual spread. Use it.

5. Halve the threshold, because p_sim is one-sided

This is the detail that trips people up, and it is easy to validate. Run the test on pure noise, where roughly 5% of units should be flagged:

noise = rng.normal(20, 5, len(zones))
gi_noise = G_Local(noise, w, star=True, permutations=999)

for label, mask in [
    ("p_sim <= 0.05", gi_noise.p_sim <= 0.05),
    ("p_sim <= 0.025", gi_noise.p_sim <= 0.025),
    ("|z_sim| > 1.96", np.abs((gi_noise.Gs - gi_noise.EG_sim) / gi_noise.seG_sim) > 1.96),
]:
    print(f"  {label:18} {mask.sum():3} of {len(noise)}  ({mask.mean():.1%})")
  p_sim <= 0.05       31 of 400  (7.8%)
  p_sim <= 0.025      19 of 400  (4.8%)
  |z_sim| > 1.96      18 of 400  (4.5%)

p_sim <= 0.05 flags well over the expected rate. p_sim <= 0.025 and |z_sim| > 1.96 both land near 5%, which is what a correctly calibrated two-sided test at Ξ± = 0.05 should do.

Run this noise control on your own weights. It takes one line and it is the only way to know your threshold is calibrated.

A noise control flagging 7.8 percent of units at p_sim under 0.05 but 4.8 percent at 0.025, matching the expected five percent.
The one-sided p-value doubles your false-positive rate unless you halve the threshold. A noise control catches it in one line.

Code examples

Example 1 β€” hotspot detection with the calibration built in 

import geopandas as gpd
import numpy as np
import pandas as pd
from esda.getisord import G_Local
from libpysal.weights import DistanceBand, KNN, Queen


def hotspots(zones, column, w, *, alpha=0.05, permutations=999, star=True):
    """Getis-Ord Gi* using permutation inference, with a two-sided threshold."""
    values = zones[column].to_numpy(dtype="float64")
    if not np.isfinite(values).all():
        raise ValueError(f"{(~np.isfinite(values)).sum()} non-finite values in {column!r}")
    if np.var(values) == 0:
        raise ValueError(f"{column!r} has zero variance β€” Gi* is undefined")
    if w.islands:
        print(f"  {len(w.islands)} island(s) have no local sum and will be excluded")

    gi = G_Local(values, w, star=star, permutations=permutations)
    z_sim = (gi.Gs - gi.EG_sim) / gi.seG_sim

    # p_sim is one-sided; halve alpha for a two-sided test
    significant = gi.p_sim <= alpha / 2

    out = zones.copy()
    out["gi_z"] = z_sim
    out["gi_p"] = gi.p_sim
    out["hotspot"] = np.where(
        ~significant, "not significant",
        np.where(z_sim > 0, "hot", "cold"),
    )
    # confidence bands, for mapping
    out["band"] = pd.cut(
        np.where(significant, z_sim, 0.0),
        bins=[-np.inf, -3.29, -2.58, -1.96, 1.96, 2.58, 3.29, np.inf],
        labels=["cold 99.9%", "cold 99%", "cold 95%", "not significant",
                "hot 95%", "hot 99%", "hot 99.9%"],
    )

    print(f"  {column}: {(out['hotspot'] == 'hot').sum()} hot, "
          f"{(out['hotspot'] == 'cold').sum()} cold, "
          f"{(out['hotspot'] == 'not significant').sum()} not significant")
    return out


def calibration_check(w, *, n_trials=1, alpha=0.05, permutations=999, seed=0):
    """Run Gi* on pure noise β€” a correct test flags about alpha of the units."""
    rng = np.random.default_rng(seed)
    rates = []
    for _ in range(n_trials):
        noise = rng.normal(0, 1, w.n)
        gi = G_Local(noise, w, star=True, permutations=permutations)
        rates.append((gi.p_sim <= alpha / 2).mean())
    print(f"  noise control: {np.mean(rates):.1%} flagged "
          f"(expected {alpha:.0%}) over {n_trials} trial(s)")
    return float(np.mean(rates))


w = Queen.from_dataframe(zones, use_index=True)
w.transform = "r"
calibration_check(w)
result = hotspots(zones, "value", w)
  noise control: 4.8% flagged (expected 5%)
  value: 101 hot, 104 cold, 195 not significant

The calibration line runs in a second and turns "I used the standard threshold" into "I verified the threshold on this weights structure".

Example 2 β€” mapping the confidence bands 

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

BANDS = ["cold 99.9%", "cold 99%", "cold 95%", "not significant",
         "hot 95%", "hot 99%", "hot 99.9%"]
COLOURS = ["#08519c", "#3182bd", "#9ecae1", "#f0f0f0",
           "#fcae91", "#fb6a4a", "#cb181d"]

fig, ax = plt.subplots(figsize=(8, 8))
result.plot(column="band", categorical=True, cmap=ListedColormap(COLOURS),
            categories=BANDS, legend=True, ax=ax,
            edgecolor="white", linewidth=0.2,
            legend_kwds={"loc": "lower right", "frameon": False})
ax.set_title(f"Getis-Ord Gi* Β· {len(result)} zones Β· queen contiguity Β· 999 permutations")
ax.set_axis_off()
fig.savefig("hotspots.png", dpi=200, bbox_inches="tight")

print(result["band"].value_counts().reindex(BANDS).to_string())
cold 99.9%          49
cold 99%            31
cold 95%            23
not significant    196
hot 95%             21
hot 99%             25
hot 99.9%           55

A diverging colour ramp centred on "not significant" is the convention, and it is the right one: the middle category is genuinely a midpoint, not a low value. Using a sequential ramp here implies that "not significant" is less than "cold", which is meaningless.

Put the weights scheme and permutation count in the title. A Gi* map without them cannot be reproduced.

Example 3 β€” how much the answer depends on the weights 

def weights_sensitivity(zones, column, schemes):
    rows = []
    for label, w in schemes.items():
        w.transform = "r"
        result = hotspots(zones, column, w)
        rows.append({
            "weights": label,
            "mean_nb": round(w.mean_neighbors, 1),
            "hot": int((result["hotspot"] == "hot").sum()),
            "cold": int((result["hotspot"] == "cold").sum()),
        })
    frame = pd.DataFrame(rows)
    print(frame.to_string(index=False))
    return frame


weights_sensitivity(zones, "value", {
    "queen": Queen.from_dataframe(zones, use_index=True),
    "knn k=4": KNN.from_dataframe(zones, k=4),
    "knn k=8": KNN.from_dataframe(zones, k=8),
    "knn k=16": KNN.from_dataframe(zones, k=16),
})
  weights  mean_nb  hot  cold
    queen      7.4  101   104
  knn k=4      4.0   76    82
  knn k=8      8.0  102   106
 knn k=16     16.0  124   138

More neighbours means a larger local sum, which is smoother and further from the global mean β€” so wider neighbourhoods find more significant units, not fewer. The hot count rises from 76 to 124 as k goes from 4 to 16.

That is the spatial scale of the analysis, chosen by you. State it, and prefer a distance band when you can name the scale in metres.

Explanation

What Gi* actually computes

For each unit i, Gi* is the sum of values in i's neighbourhood (including i itself) divided by the sum of all values, compared against what that ratio would be if the values were distributed at random.

A high Gi* means the neighbourhood's total is disproportionately large. Crucially it is about the local sum, not the local value β€” a moderate unit surrounded by high ones is part of a hotspot, and a very high isolated unit may not be.

That difference is why Gi* and a plain choropleth of the values look different, and why Gi* is the better answer to "where is the problem concentrated".

Why the analytical variance fails so badly

The published variance formula for Gi* is derived under randomisation with an assumption that the local sums are approximately normal. Two things break it in practice:

  • Skewed values. Counts and rates are rarely normal; a few large values dominate.
  • Small neighbourhoods. A sum of eight values is not close to normal regardless of the underlying distribution.

Both inflate the formula's variance estimate relative to the true sampling distribution. On the data above the ratio was 3.11, which flattened the z-scores from a genuine range of Β±5.2 to a useless Β±1.8.

Conditional permutation sidesteps the whole issue by measuring the distribution empirically. It costs 999 shuffles β€” a fraction of a second β€” and it is the reason p_sim exists.

Why the one-sided p-value needs halving

esda computes p_sim as the proportion of permuted values at least as extreme as the observed one in the observed direction, which caps at 0.5. It is a one-tailed p-value.

Hotspot analysis is two-tailed: you care about both unusually high and unusually low neighbourhoods. Comparing a one-tailed p against 0.05 therefore tests at an effective Ξ± of 0.10, and the noise control shows exactly that β€” 7.8% flagged instead of 5%.

Halving the threshold to 0.025 restores calibration (4.8% on noise). The equivalent test on the permutation z-score, |z_sim| > 1.96, gives 4.5% β€” the same answer by a different route.

A moderate unit surrounded by high neighbours forming part of a hotspot, and an isolated high unit that does not.
Gi* measures the neighbourhood total, not the unit's own value. That is what separates it from a choropleth.

Why multiple comparisons still matter

Even with a calibrated threshold you are running one test per unit. At 400 units and Ξ± = 0.05 you expect about 20 false positives, and the noise control found 19.

For a map, that is usually acceptable β€” false positives are scattered, while genuine hotspots are contiguous blocks, and the eye discounts isolated cells. For a claim about a specific unit, it is not. Apply a false-discovery-rate correction, and note that the tests are not independent (neighbouring units share data), so standard corrections are conservative in a way that is hard to quantify.

The pragmatic position: report the threshold, report the noise-control rate, and treat isolated single-cell hotspots with suspicion.

Edge cases or notes

  • Never use Zs or p_norm. On the example data they found zero hotspots where permutation inference found 205.
  • star=True includes the unit itself. star=False answers a different question and is rarely what "hotspot" means.
  • Permutation results move between runs. Set np.random.seed() before G_Local if the exact counts must be reproducible; the split moves by a few units otherwise.
  • 999 permutations gives a minimum p_sim of 0.001. For stricter thresholds after correction, use 9,999.
  • Gi* needs a variable with meaningful magnitude β€” counts, rates, amounts. It is not appropriate for categorical or purely ordinal data.
  • Islands have no local sum. With a distance band this is common; check w.islands before running.
  • Row-standardise for irregular units so a unit with twelve neighbours is not weighted three times a unit with four.
  • Hotspots move with the zone size. A Gi* map is a statement about one aggregation β€” see the modifiable areal unit problem.
  • Local Moran (LISA) answers a different question: it finds both clusters and spatial outliers, where Gi* only finds high and low concentrations.

FAQ

Why did Gi* find no hotspots in obviously clustered data?

You used Zs. Its analytical variance is unreliable β€” over three times too large on the example here β€” so nothing reaches 1.96. Use p_sim or the permutation z-score instead.

What threshold should I use with p_sim?

Half your intended Ξ±, because p_sim is one-sided. For a two-sided test at 5%, use p_sim <= 0.025. Verify with a noise control.

What does star=True do?

Includes the unit itself in its own neighbourhood, making it Gi* rather than Gi. It is what you want in almost every case.

How do I know my threshold is calibrated?

Run Gi* on random noise with the same weights. A correct 5% test flags about 5% of units. If it flags 9%, your threshold is one-sided.

Which weights should I use?

A distance band when you can name the spatial scale in metres β€” it makes the scale explicit. Contiguity or KNN otherwise. Wider neighbourhoods find more significant units, so report the choice.

What is the difference between Gi* and local Moran's I?

Gi* finds concentrations of high or low values. Local Moran also identifies spatial outliers β€” a low unit among high neighbours β€” which Gi* cannot express.

Do I need to correct for multiple comparisons?

For claims about specific units, yes. For a map, the convention is to report the threshold and the noise-control rate, and to treat isolated single-cell hotspots as suspect.