Run a Hot Spot Analysis with Getis-Ord Gi* in PyQGIS

A density map shows where events concentrate, but every random scatter of points has some dense areas. A hot spot analysis asks a sharper question: is the concentration around this location higher than you would expect if values were spread at random across the study area? The Getis-Ord Gi* statistic answers it for every feature, returning a z-score — large and positive for hot spots, large and negative for cold spots, near zero for nothing unusual.

This recipe belongs to Spatial Statistics & Pattern Analysis. QGIS has no built-in Gi* algorithm, so the recipe computes it directly: aggregate points into a regular grid, define which cells are neighbours, calculate Gi* with NumPy, correct for testing many cells at once, and map hot and cold spots with a diverging style.

From points to significant hot spotsIncident points are counted into hexagon cells. Each cell's neighbourhood is defined, for example all cells within 600 metres including itself. Gi* compares the sum of counts in each neighbourhood with what the overall mean would predict, scaled by the variance, giving a z-score. Cells with z-scores beyond the corrected significance threshold are hot or cold spots.Count, define neighbours, score, testpointsincidentshexagon countsone value per cellneighbourhoodsdistance bandincludes selfGi* z-scoreper cellsignificant after correction → hot · cold · not significant

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series. NumPy ships with QGIS.
  • A point layer of events in a projected CRS, and the boundary of the study area.
  • A neighbourhood scale in mind — the distance at which you expect events to influence each other.

Aggregate points to a grid

Gi* works on values attached to areas. Counting points into a regular hexagon grid gives every cell a value and every cell the same size and shape, which keeps neighbourhoods comparable.

import processing
from qgis.core import QgsProject

incidents = QgsProject.instance().mapLayersByName("burglaries")[0]
city = QgsProject.instance().mapLayersByName("city_boundary")[0]

grid = processing.run("native:creategrid", {
    "TYPE": 4, "EXTENT": city.extent(), "HSPACING": 250, "VSPACING": 250,
    "CRS": city.crs(), "OUTPUT": "memory:"})["OUTPUT"]
grid = processing.run("native:extractbylocation", {
    "INPUT": grid, "PREDICATE": [0], "INTERSECT": city, "OUTPUT": "memory:"})["OUTPUT"]
counts = processing.run("native:countpointsinpolygon", {
    "POLYGONS": grid, "POINTS": incidents, "FIELD": "n",
    "OUTPUT": "memory:hex_counts"})["OUTPUT"]
print(counts.featureCount(), "cells;", sum(f["n"] for f in counts.getFeatures()), "points")

Breakdown: Grid type 4 is hexagonal; hexagons have six equidistant neighbours, which suits distance-based neighbourhoods better than squares. Keeping only cells that intersect the city removes empty sea or countryside cells that would drag the mean down and make everything look hot. Cell size should be small enough to locate hot spots usefully but large enough that most cells contain a few events; creating a hexagon grid and counting points covers the choice.

Define neighbours

Gi* sums values over each cell's neighbourhood, so the neighbourhood definition is the analysis's scale. A distance band — all cells whose centres lie within a given distance, including the cell itself — is the standard choice.

A distance band neighbourhoodA central hexagon and the ring of cells whose centres lie within the distance band, here 600 metres, around it. All cells inside the band, including the central one, form the neighbourhood with weight one; cells outside have weight zero. A larger band smooths results and finds broader hot spots; a smaller band finds tighter ones.Everything within the band counts, including itselfthe cell itselfweight 1 (the * in Gi*)centres within 600 mweight 1outside the bandweight 0

import numpy as np
from qgis.core import QgsSpatialIndex

feats = list(counts.getFeatures())
ids = [f.id() for f in feats]
x = np.array([f["n"] for f in feats], dtype=float)
centres = [f.geometry().centroid() for f in feats]
pos = {fid: i for i, fid in enumerate(ids)}

BAND = 600  # metres
index = QgsSpatialIndex()
for f, c in zip(feats, centres):
    index.addFeature(f.id(), c.boundingBox())

neighbours = []
for i, c in enumerate(centres):
    cand = index.intersects(c.boundingBox().buffered(BAND))
    neighbours.append([pos[j] for j in cand if c.distance(centres[pos[j]]) <= BAND])
sizes = [len(n) for n in neighbours]
print(f"neighbours per cell: min {min(sizes)}, median {sorted(sizes)[len(sizes) // 2]}, max {max(sizes)}")

Breakdown: Indexing cell centres rather than hexagons makes the distance test a centre-to-centre distance, which is the conventional definition. Each list includes the cell itself, because its distance to itself is zero — that is what distinguishes Gi* from Gi. Every cell should have several neighbours; cells with only themselves give meaningless z-scores, so if the minimum is 1, the band is too small for the grid spacing. A band of two to three cell widths is a reasonable start, and comparing results at two bands shows how scale-dependent the hot spots are.

Compute Gi* with NumPy

With values and neighbour lists in hand, the statistic is a few lines of arithmetic: compare each neighbourhood's sum with what the global mean predicts, and scale by the expected spread.

n = len(x)
xbar = x.mean()
s = np.sqrt((x ** 2).mean() - xbar ** 2)

z = np.zeros(n)
for i, nb in enumerate(neighbours):
    w_sum = len(nb)                         # binary weights: Σw = Σw²
    local = x[nb].sum()
    num = local - xbar * w_sum
    den = s * np.sqrt((n * w_sum - w_sum ** 2) / (n - 1))
    z[i] = num / den if den > 0 else 0.0
print(f"z range {z.min():.2f} to {z.max():.2f}")

Breakdown: With binary weights — 1 inside the band, 0 outside — the sum of weights and the sum of squared weights are both just the neighbour count, which simplifies the formula. The numerator is how much the neighbourhood's total exceeds what the global mean would give for that many cells; the denominator is its standard deviation under randomness. The result is a z-score directly: values beyond ±1.96 are significant at the 5 % level for a single test — but this analysis runs hundreds of tests at once, which the next step addresses. If s is zero every cell has the same count and there is nothing to find.

Correct for many tests

Testing every cell at the 5 % level guarantees false hot spots: with a thousand cells, about fifty pass by chance alone. Neighbouring cells also share data, so their tests are not independent. A false discovery rate correction is the usual remedy.

Uncorrected versus corrected thresholdsWith one thousand cells tested at the five percent level, roughly fifty appear significant by chance alone. The Benjamini–Hochberg false discovery rate correction ranks p-values and raises the bar so that, among cells called significant, the expected share of false discoveries stays at five percent. Fewer cells pass, and those that do are trustworthy.Many tests need a stricter baruncorrected, |z| > 1.961000 cells × 5 % ≈ 50chance hot spotsFDR correctedexpected false share 5 %fewer, trustworthy

from math import erf, sqrt

p = np.array([2 * (1 - 0.5 * (1 + erf(abs(v) / sqrt(2)))) for v in z])   # two-sided
order = np.argsort(p)
m = len(p)
thresh = 0.05 * (np.arange(1, m + 1) / m)
passed = p[order] <= thresh
cutoff = p[order][passed].max() if passed.any() else 0.0
significant = p <= cutoff

hot = int(((z > 0) & significant).sum())
cold = int(((z < 0) & significant).sum())
print(f"{hot} hot and {cold} cold cells after FDR correction (cutoff p = {cutoff:.4f})")

Breakdown: The two-sided p-value comes from the standard normal distribution via the error function, so no SciPy is needed. The Benjamini–Hochberg procedure sorts p-values, compares each with an increasing threshold, and declares significant everything up to the largest p-value that passes. It controls the expected share of false discoveries among the cells called significant, which is the property a hot spot map needs. Bonferroni — dividing 0.05 by the number of tests — is simpler but far stricter and usually finds nothing on real data.

Write results and map them

Writing z-scores, p-values and a class back to the grid makes the result a normal layer: style it, label it, export it.

from qgis.core import QgsField, QgsCategorizedSymbolRenderer, QgsRendererCategory, QgsFillSymbol
from qgis.PyQt.QtCore import QVariant

prov = counts.dataProvider()
prov.addAttributes([QgsField("gi_z", QVariant.Double), QgsField("gi_p", QVariant.Double),
                    QgsField("spot", QVariant.String, len=8)])
counts.updateFields()
iz, ip, isp = (counts.fields().indexOf(n) for n in ("gi_z", "gi_p", "spot"))
changes = {}
for i, fid in enumerate(ids):
    label = "hot" if significant[i] and z[i] > 0 else "cold" if significant[i] else "none"
    changes[fid] = {iz: float(z[i]), ip: float(p[i]), isp: label}
prov.changeAttributeValues(changes)

cats = [QgsRendererCategory(v, QgsFillSymbol.createSimple({"color": c, "outline_style": "no"}), v)
        for v, c in (("hot", "#b91c1c"), ("cold", "#2563eb"), ("none", "#e8e4d8"))]
counts.setRenderer(QgsCategorizedSymbolRenderer("spot", cats))
QgsProject.instance().addMapLayer(counts)

Breakdown: Bulk changeAttributeValues writes every cell in one provider call. Three categories — hot, cold, not significant — make the map honest: a graduated ramp on raw z-scores invites readers to see hot spots in cells that did not pass the test. Red for hot and blue for cold is the convention; a neutral fill for the rest keeps attention on what matters. For the z-scores themselves, a diverging ramp centred on zero is the right choice if you show them at all — see graduated renderers.

Check sensitivity before you report

Hot spots that appear at one cell size and band but vanish at another are fragile. A quick sensitivity check reruns the analysis over a few settings and keeps the cells that are hot in most of them.

def gi_star(values, centres, band, q=0.05):
    n = len(values)
    xbar, sd = values.mean(), values.std()
    zs = np.zeros(n)
    for i, c in enumerate(centres):
        cand = index.intersects(c.boundingBox().buffered(band))
        nb = [pos[j] for j in cand if c.distance(centres[pos[j]]) <= band]
        w = len(nb)
        den = sd * np.sqrt((n * w - w * w) / (n - 1))
        zs[i] = (values[nb].sum() - xbar * w) / den if den > 0 else 0.0
    ps = np.array([2 * (1 - 0.5 * (1 + erf(abs(v) / sqrt(2)))) for v in zs])
    srt = np.sort(ps)
    ok = srt <= q * np.arange(1, n + 1) / n
    cut = srt[ok].max() if ok.any() else 0.0
    return zs, ps <= cut

stable = np.zeros(n, dtype=int)
for band in (450, 600, 900):
    zb, sig = gi_star(x, centres, band)
    stable += ((zb > 0) & sig).astype(int)
print((stable >= 2).sum(), "cells hot at two or more bands")

Breakdown: Wrapping the neighbour construction, z-scores and correction in one function makes reruns trivial; values.std() is the same population standard deviation as s above, and the spatial index of cell centres is reused across bands. Cells hot at most bands are robust findings worth acting on; cells hot at one band only are scale artefacts. Report the band and cell size with every hot spot map, because a reader cannot judge the result without them.

QGIS version compatibility

The algorithms used — native:creategrid, native:extractbylocation, native:countpointsinpolygon — have the same parameters on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. On QGIS 4, QgsField takes QMetaType.Type.Double and QMetaType.Type.QString instead of QVariant types. NumPy is bundled with all official QGIS installers.

Troubleshooting

  • Everything near the edge is cold. Empty cells outside the real area were kept; clip the grid to the study area.
  • No cells are significant after correction. The band is too small, the grid too fine, or there genuinely is no hot spot.
  • z-scores are all NaN. Every cell has the same value, so the standard deviation is zero.
  • Hot spots follow population. Counts reflect where people live; analyse rates (events per household) instead.

Conclusion

Count events into a hexagon grid clipped to the study area, define neighbours with a distance band that includes each cell itself, compute Gi* z-scores with NumPy, correct for many tests with the false discovery rate, map hot, cold and not significant as three categories, and report only spots that survive a change of scale.

Frequently Asked Questions

Can I use PySAL instead? Yes. The esda package's G_Local with star=True computes the same statistic with permutation-based p-values; install it into QGIS's Python first.

Should I analyse counts or rates? Rates, when the population at risk varies — otherwise hot spots simply show where people are.

What about Local Moran's I? It finds clusters of similar values and outliers (high surrounded by low); Gi* only finds concentrations of high or low values. Use Moran's I when outliers matter.

Does the cell shape matter? Hexagons reduce edge-orientation effects; squares work but have uneven neighbour distances.