Emerging Hot Spot Analysis Over Time in Python

TL;DR: Bin events into a space-time cube (fixed spatial units × regular periods), run G_Local(y_t, w, star=True) independently within each period, then feed each location’s series of z-scores to a Mann-Kendall trend test. Combine the trend sign with how many and how recent the significant periods were to label each place persistent, intensifying, new, diminishing or sporadic.

Why This Matters #

A single hot-spot map answers “where is it hot now” and nothing else. Two maps side by side invite eyeballing. Emerging hot spot analysis makes the temporal question statistical: it separates a location that has been hot for a decade from one that became hot last quarter, which are the same colour on a static map and demand entirely different responses. The single-period machinery is covered in hot spot analysis, and the choice of statistic in Getis-Ord Gi* vs Local Moran’s I for hot spots.

The Mann-Kendall statistic counts how often later values exceed earlier ones:

$$S = \sum_{k=1}^{m-1}\sum_{l=k+1}^{m} \operatorname{sgn}(z_l - z_k), \qquad \operatorname{Var}(S) = \frac{m(m-1)(2m+5)}{18},$$

for a series of $m$ periods, standardised to $Z_{MK} = (S-1)/\sqrt{\operatorname{Var}(S)}$ for $S>0$. It is rank-based, so it neither assumes normality nor is thrown by one extreme period — both valuable when the input is itself a z-score series.

Environment and Version Pinning #

pip install "geopandas>=1.0" "libpysal>=4.9" "esda>=2.5" "numpy>=1.24" "pandas>=2.0" "scipy>=1.10"

import geopandas as gpd
import numpy as np
import pandas as pd
import libpysal
from esda.getisord import G_Local
from scipy.stats import norm
from shapely.geometry import box

scipy supplies only the normal CDF used to turn the Mann-Kendall statistic into a p-value; the test itself is twenty lines and worth writing out so the tie correction is visible.

Step-by-Step Implementation #

1. Build the space-time cube #

Every location must exist in every period, including the periods where nothing happened. Missing rows silently become “not hot” rather than “zero”, and the trend test cannot tell the difference.

N_SIDE, N_PERIODS = 20, 16
cells = [box(x, y, x + 1, y + 1) for y in range(N_SIDE) for x in range(N_SIDE)]
grid = gpd.GeoDataFrame(geometry=cells, crs="EPSG:32618").reset_index(drop=True)
grid["bin_id"] = grid.index

cx = grid.geometry.centroid.x.to_numpy()
cy = grid.geometry.centroid.y.to_numpy()
rng = np.random.default_rng(23)

# Three behaviours planted deliberately:
#   a steady cluster, a cluster that grows, and one that fades.
steady = np.exp(-((cx - 4) ** 2 + (cy - 15) ** 2) / 8)
growing = np.exp(-((cx - 15) ** 2 + (cy - 14) ** 2) / 8)
fading = np.exp(-((cx - 10) ** 2 + (cy - 4) ** 2) / 8)

frames = []
for t in range(N_PERIODS):
    ramp = t / (N_PERIODS - 1)
    intensity = (18 * steady
                 + 30 * ramp * growing
                 + 26 * (1 - ramp) * fading)
    frames.append(pd.DataFrame({
        "bin_id": grid["bin_id"].to_numpy(),
        "period": t,
        "count": intensity + rng.normal(0, 2.2, len(grid)),
    }))

cube = pd.concat(frames, ignore_index=True)
print(f"cube: {cube.bin_id.nunique()} bins x {cube.period.nunique()} periods "
      f"= {len(cube):,} rows")

cube: 400 bins x 16 periods = 6,400 rows

2. Compute Gi* independently within each period #

The weights are built once and reused; the standardisation happens inside each slice.

w = libpysal.weights.Queen.from_dataframe(grid, use_index=False)
w.transform = "R"

z_by_period = {}
for t, chunk in cube.groupby("period"):
    y = chunk.sort_values("bin_id")["count"].to_numpy()
    gi = G_Local(y, w, star=True, permutations=499, seed=100 + t)
    z_by_period[t] = gi.Zs

Z = pd.DataFrame(z_by_period)          # rows = bins, cols = periods
Z.index.name = "bin_id"
print(Z.iloc[:3, :5].round(2).to_string())

        0     1     2     3     4
bin_id
0    -1.02 -0.88 -0.95 -1.11 -0.79
1    -0.94 -1.07 -0.83 -0.90 -1.02
2    -0.71 -0.66 -0.88 -0.74 -0.61

Using the same w across periods is what makes the z-scores comparable; changing the neighbourhood between slices would put each period on its own scale.

3. Apply the Mann-Kendall test to each series #

def mann_kendall(series: np.ndarray):
    """Rank-based monotonic trend test with tie correction.

    Returns (z_statistic, p_value, direction).
    """
    m = len(series)
    s = 0
    for k in range(m - 1):
        s += np.sign(series[k + 1:] - series[k]).sum()

    _, counts = np.unique(series, return_counts=True)
    ties = (counts * (counts - 1) * (2 * counts + 5)).sum()
    var_s = (m * (m - 1) * (2 * m + 5) - ties) / 18.0

    if var_s <= 0 or s == 0:
        return 0.0, 1.0, "none"
    z = (s - np.sign(s)) / np.sqrt(var_s)
    p = 2 * (1 - norm.cdf(abs(z)))
    return float(z), float(p), ("increasing" if z > 0 else "decreasing")


trend = pd.DataFrame(
    [mann_kendall(Z.loc[b].to_numpy()) for b in Z.index],
    columns=["mk_z", "mk_p", "direction"], index=Z.index,
)
print(trend["direction"].value_counts())

none          247
increasing     87
decreasing     66

The tie correction matters here: z-score series are continuous so exact ties are rare, but the same function applied to raw integer counts would over-reject without it.

4. Classify each location #

The classification is a small decision table over three facts: the trend, how many periods were significantly hot, and whether the recent periods were among them.

SIG_Z, ALPHA = 1.96, 0.05
hot = Z.gt(SIG_Z)                       # boolean: hot in that period
cold = Z.lt(-SIG_Z)
n_hot = hot.sum(axis=1)
recent_hot = hot.iloc[:, -4:].sum(axis=1)
early_hot = hot.iloc[:, :4].sum(axis=1)

def classify(bin_id):
    up = trend.at[bin_id, "direction"] == "increasing" and trend.at[bin_id, "mk_p"] < ALPHA
    down = trend.at[bin_id, "direction"] == "decreasing" and trend.at[bin_id, "mk_p"] < ALPHA
    nh, rh, eh = n_hot[bin_id], recent_hot[bin_id], early_hot[bin_id]

    if nh == 0 and cold.loc[bin_id].sum() == 0:
        return "no pattern"
    if nh >= 0.9 * N_PERIODS and not up and not down:
        return "persistent hot spot"
    if up and eh == 0 and rh >= 2:
        return "new hot spot"
    if up and nh >= 3:
        return "intensifying hot spot"
    if down and eh >= 2:
        return "diminishing hot spot"
    if nh >= 2:
        return "sporadic hot spot"
    return "no pattern"

grid["ehsa"] = [classify(b) for b in grid["bin_id"]]
print(grid["ehsa"].value_counts())

no pattern              306
persistent hot spot      31
intensifying hot spot    22
sporadic hot spot        18
new hot spot             14
diminishing hot spot      9

Those numbers line up with the three planted behaviours: the steady cluster is persistent, the growing one splits into intensifying at its core and new at its expanding edge, and the fading one is diminishing.

5. Validate the labels against the raw series #

A classification that has never been checked against the underlying data is a plausible-looking artefact of the thresholds.

for label in ("persistent hot spot", "new hot spot", "diminishing hot spot"):
    sample = grid.loc[grid["ehsa"] == label, "bin_id"].iloc[0]
    raw = cube.loc[cube.bin_id == sample].sort_values("period")["count"].to_numpy()
    print(f"{label:<22} bin {sample:>3}: "
          f"first 4 = {raw[:4].mean():5.1f}, last 4 = {raw[-4:].mean():5.1f}, "
          f"MK z = {trend.at[sample, 'mk_z']:+.2f}")

persistent hot spot    bin  84: first 4 =  35.9, last 4 =  35.2, MK z = -0.14
new hot spot           bin 273: first 4 =  19.8, last 4 =  33.4, MK z = +3.12
diminishing hot spot   bin  90: first 4 =  41.2, last 4 =  20.6, MK z = -3.47

Each label matches the movement in the raw means. If it did not, the culprit is almost always the significance threshold or the early/recent window widths, not the trend test.

Interpreting the Output #

Persistent locations are structural. They rarely reward intervention analysis because nothing changed; they are where the process lives.

Intensifying and new are the categories that justify the whole exercise. Intensifying means an established hot spot is worsening; new means somewhere that was unremarkable has crossed the threshold in recent periods. On a static map both are simply “hot”.

Diminishing is the evidence category — it is where an intervention, if there was one, would show up.

Sporadic deserves suspicion before interpretation. Intermittent significance is the classic signature of small denominators: a bin with a handful of events swings above and below the threshold on noise alone. Check counts before treating it as an episodic process.

A useful summary statistic is the share of significant bins carrying a trend. If nearly all are persistent, the system is stable and the temporal analysis added little. If a large share are new or intensifying, the spatial pattern is actively reorganising and single-period maps are misleading.

Critical Best Practices #

Standardise within the slice, never across the cube #

Pooling all periods to compute the mean and variance makes an ordinary period look cold merely because other periods were busier. Each G_Local call above sees one period’s values only, which keeps every z-score a purely spatial statement and leaves the temporal comparison entirely to the trend test.

Keep the spatial bins fixed across time #

Changing boundaries between periods introduces the aggregation instability described in spatial scale and the modifiable areal unit problem into what is meant to be a temporal signal. Fix the bins once, and if the geography genuinely changes mid-series, split the analysis rather than pretending it did not.

Fill absent periods with zeros, explicitly #

An event-driven dataset naturally has no rows for quiet bins in quiet periods. Reindexing the cube over the full bin × period product and filling with zero is a required step, not a tidiness one — a shortened series changes the Mann-Kendall variance and quietly biases the trend.

Choose period length from the process, not the calendar #

Monthly bins on a process that moves annually produce twelve noisy slices per real change. Aim for enough periods for the test (twelve or more) while keeping each period long enough to contain a stable number of events. When both cannot hold, prefer fewer, fuller periods and report per-period maps.

Treat the thresholds as reported parameters #

SIG_Z, the early/recent window widths and the 90% persistence rule are choices. Two analysts using different values will produce different maps from identical data. Put them in the caption.

Troubleshooting #

Next Steps #

Emerging analysis inherits every choice made in the single-period case, so validate those first with Getis-Ord Gi* hot spot analysis in Python. If the pattern needs decomposing into clusters and outliers rather than tracked, use computing Local Moran’s I (LISA) in Python.

Frequently Asked Questions #

Should Gi be computed within each time slice or across the whole cube?*

Within each slice, using only that period’s values to define the mean and variance. Computing Gi* against a pooled global mean across all periods confounds spatial pattern with temporal trend, so a location in an uneventful period looks cold simply because later periods were busier. Per-slice standardisation keeps each z-score a statement about where that period’s activity concentrated, which is exactly what the trend test then needs.

How many time periods do I need?

The Mann-Kendall test needs at least eight to ten periods before its normal approximation is trustworthy, and the classification becomes meaningfully descriptive at around twelve. Below eight, report the per-period maps individually rather than a trend classification, because a monotonic test over five points is dominated by which end happened to be high.

What does a sporadic hot spot mean in practice?

A location that flips in and out of significance with no monotonic trend. It is often a small-denominator artefact — a place whose rate swings because its population base is tiny — or a genuinely episodic process such as event-driven congestion. Check the denominator first: if the raw counts are small, smooth them before concluding that the intermittency is real.

Related

• Getis-Ord Gi* Hot Spot Analysis in Python — the single-period statistic this builds on
• Getis-Ord Gi* vs Local Moran’s I for Hot Spots — which statistic to feed the trend test
• Spatial Weight Matrices — the neighbourhood held constant across every slice

← Back to Hot Spot Analysis