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
from libpysal.weights import Queen

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

10 hotspots

Ten, on the crime rates of Leeds' 488 neighbourhoods (see The data used on this page), which have a Moran's I of 0.35 at p = 0.001. And no cold spots at all: the analytical z-scores run from −0.81 to +5.11, so nothing can ever fall below −1.96.

The problem is Zs. On this data its analytical standard deviation is 2.8 times larger than the permutation-based one. The permutation test finds 19 hot spots and 31 cold ones.

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["rate"].to_numpy(), w, star=True, permutations=999, seed=0)

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
not significant    438
cold                31
hot                 19

Areal units, weights, star=True, permutation p-values

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

The examples here use the Leeds LSOAs (see The data used on this page), so hexes below is simply the zones:

hexes = zones

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=488, mean neighbours 5.7, 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.

band = DistanceBand.from_dataframe(hexes, threshold=800, binary=True, silence_warnings=True)
print(f"800 m band: mean neighbours {band.mean_neighbors:.1f}, islands {len(band.islands)} of {band.n}")

800 m band: mean neighbours 3.3, islands 79 of 488

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

values = hexes["rate"].to_numpy(dtype="float64")
gi = G_Local(values, w, star=True, permutations=999, seed=0)

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

z_sim = (gi.Gs - gi.EG_sim) / gi.seG_sim
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  -0.81 .. +5.11
permutation z range  -1.77 .. +12.14
analytic sd / sim sd 2.83

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:

rng = np.random.default_rng(0)
noise = rng.normal(20, 5, len(zones))
gi_noise = G_Local(noise, w, star=True, permutations=999, seed=0)

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       39 of 488  (8.0%)
  p_sim <= 0.025      25 of 488  (5.1%)
  |z_sim| > 1.96      26 of 488  (5.3%)

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.

Calibrated hotspots and confidence bands for Leeds

The data used on this page

Every example runs on real data: the 488 Lower Layer Super Output Areas (LSOAs) of Leeds from the ONS, their Census 2021 population from Nomis, and the street-level crimes West Yorkshire Police published for June 2026 on data.police.uk. rate is crimes per 1,000 residents.

import io
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query")
features, offset = [], 0
while True:                                   # the service pages at 1,000 features
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD,LSOA21NM",
                                     "outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700).set_index("LSOA21CD").sort_index()

population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(zones.index)}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = population

west, south, east, north = zones.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                  data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if reply.status_code != 429:
                break
            time.sleep(2 ** attempt)
        reply.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "lon": float(c["location"]["longitude"]),
                    "lat": float(c["location"]["latitude"])} for c in reply.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes.lon, crimes.lat), crs=4326).to_crs(27700)
incidents = gpd.sjoin(incidents, zones[["geometry"]], predicate="within").rename(columns={"index_right": "lsoa", "LSOA21CD": "lsoa"})
zones["crimes"] = incidents["lsoa"].value_counts().reindex(zones.index, fill_value=0)
zones["rate"] = 1000 * zones["crimes"] / zones["population"]
print(f"{len(zones)} LSOAs, {zones['population'].sum():,} residents, {len(incidents):,} crimes; "
      f"rate median {zones['rate'].median():.1f}, max {zones['rate'].max():.1f} per 1,000")

488 LSOAs, 811,957 residents, 9,229 crimes; rate median 8.4, max 174.4 per 1,000

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, seed=0):
    """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, seed=seed)   # seeded: counts repeat
    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 — from the permutation p-value, not from z thresholds:
    # on skewed data a significant cold unit can have z_sim above -1.96
    level = np.select([gi.p_sim <= 0.001, gi.p_sim <= 0.005, gi.p_sim <= alpha / 2], ["99.9%", "99%", "95%"], "")
    out["band"] = np.where(level == "", "not significant", np.where(z_sim > 0, "hot ", "cold ") + level)

    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, seed=seed)
        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, "rate", w)

  noise control: 5.1% flagged (expected 5%) over 1 trial(s)
  rate: 19 hot, 31 cold, 438 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())

band
cold 99.9%           2
cold 99%             9
cold 95%            20
not significant    438
hot 95%             10
hot 99%              2
hot 99.9%            7

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, "rate", {
    "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),
})

  rate: 19 hot, 31 cold, 438 not significant
  rate: 19 hot, 30 cold, 439 not significant
  rate: 14 hot, 53 cold, 421 not significant
  rate: 26 hot, 96 cold, 366 not significant
 weights  mean_nb  hot  cold
   queen      5.7   19    31
 knn k=4      4.0   19    30
 knn k=8      8.0   14    53
knn k=16     16.0   26    96

More neighbours means a larger local sum, which is smoother and more stable — so wider neighbourhoods find more significant units, not fewer. The cold count rises from 30 to 96 as k goes from 4 to 16; the hot count stays between 14 and 26, because Leeds' high rates sit in a few small central areas.

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.

What Gi* computes and why the analytical variance fails

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 2.83. Crime rates are also skewed — a few city-centre LSOAs have rates above 100 per 1,000 — so the analytical z-scores were squeezed on the low side, from −1.77 to −0.81, and no unit could ever be called cold.

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 — 8.0% flagged instead of 5%.

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

Why multiple comparisons still matter

Even with a calibrated threshold you are running one test per unit. At 488 units and α = 0.05 you expect about 24 false positives, and the noise control found 25.

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.

Zs, p_norm, islands and moving hotspots

• Never use Zs or p_norm. On the Leeds data they found 10 hot spots and could not find a cold one; permutation inference found 19 and 31.
• 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.

Internal links

• Spatial autocorrelation explained — the global statistic, and the LISA alternative
• How to calculate Moran's I in Python — establishing there is structure before locating it
• Spatial weights warn about islands — the weights failures that break Gi*
• Kernel density explained — a description where Gi* is a test
• How to bin points into hexagons in Python — producing the units Gi* needs
• The modifiable areal unit problem explained — why the zone size changes the hotspots
• Choropleth classification explained — why the diverging ramp is not optional here
• How to make a point density heatmap in Python — the continuous view of the same data
• How to map forest loss hotspots with H3 — Gi* on a hexagon grid of deforestation

FAQ

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

You used Zs. Its analytical variance is unreliable — 2.8 times too large on the Leeds example — so few units reach 1.96, and on skewed data the low side may never reach −1.96 at all. 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.