Optimizing Bandwidth for Getis-Ord Gi* Heatmaps

This guide solves one narrow problem: choosing the spatial weights bandwidth for a Getis-Ord Gi* heatmap by metric-driven search instead of visual guesswork, so the hotspots a response team maps are statistically defensible rather than an artifact of an arbitrary distance threshold. It is part of the Getis-Ord Gi* Hotspot Detection method family within the broader Disease Clustering & Spatial Statistical Modeling pipeline.

Problem Context & Constraints

Most Gi* tutorials hardcode a single distance band — DistanceBand.from_dataframe(gdf, threshold=5000) — and plot whatever falls out. That default path breaks the moment surveillance data hits it, for four reasons specific to public-health hotspot mapping.

Bandwidth is the spatial weights matrix. Gi* compares a local weighted sum of values against the global expectation under spatial randomness; the weights that define “local” come entirely from the bandwidth. A 2 km threshold and an 8 km threshold computed on the same incidence layer produce different neighbor sets, different z-scores, and different significant features. The bandwidth is the single most consequential parameter on the page, yet it is the one most often picked by eye.

Heterogeneous population density poisons a fixed threshold. A distance band that gives a dense urban tract 40 neighbors leaves a rural tract with zero. Over-connected urban cores regress toward the global mean and wash out real micro-clusters, while isolated rural features fall below the minimum-neighbor floor and either error out or return meaningless statistics. A single global distance cannot serve both regimes.

The naive map has no significance discipline. Plotting raw Gi* z-scores invites eyeballing “hot” colors as findings. With hundreds or thousands of simultaneous local tests, uncorrected p-values guarantee false hotspots. Bandwidth selection must be judged on the count of features that survive False Discovery Rate correction, not on raw color intensity.

The response curve is single-peaked, so the optimum is findable. As bandwidth widens, FDR-significant hotspot yield rises out of an under-connected, noisy regime, peaks where local structure is sharpest, then falls as wide neighborhoods dilute local sums toward the global mean and clusters dissolve. Because the curve has one interior maximum, a bounded sweep over the metric distance range will locate it — provided the search is scored on stabilized significance rather than visual preference.

Prerequisites

This procedure assumes a specific, pinned stack and an analysis-ready input layer. Mismatched versions or an un-projected CRS are the two most common causes of silently wrong bandwidths.

• Library versions (pin these in requirements.txt): geopandas==0.14.4, libpysal==4.12.1, esda==2.5.1, statsmodels==0.14.2, numpy==1.26.4, pandas==2.2.2. The esda API used below — G_Local(..., star=True), .Zs, .p_sim — is stable across the 2.5.x line.
• CRS state: the input must be in a projected, metre-based CRS before any distance threshold is computed. A bandwidth of 5000 is meaningless in EPSG:4326 degrees. Enforce the projection step described in Coordinate Reference Systems for Public Health — the code below auto-projects to the estimated UTM zone as a guard, but explicit EPSG selection is preferred for production.
• Input data: a GeoDataFrame of stable areal units (census tracts, ZCTAs, block groups) carrying a rate-based analysis variable — an age-standardized incidence rate or standardized morbidity ratio, never raw case counts. Direct patient coordinates must already have been aggregated and de-identified upstream.
• No NaN/Inf in the analysis variable. Impute or filter before the sweep; the function raises on contaminated input rather than propagating garbage z-scores.

Step-by-Step Solution

The sweep below walks a bounded metric distance range, builds a row-standardized distance-band weights matrix at each step, computes Gi*, applies Benjamini-Hochberg FDR, scores each candidate on stabilized hotspot yield, and returns the winning bandwidth plus an audit-ready optimization log. The same weights-construction and FDR conventions are shared with Getis-Ord Gi* Hotspot Detection, so a bandwidth chosen here drops straight into the production scoring pass.

# Pinned: geopandas==0.14.4, libpysal==4.12.1, esda==2.5.1,
#         statsmodels==0.14.2, numpy==1.26.4, pandas==2.2.2
import geopandas as gpd
import libpysal
import numpy as np
import pandas as pd
from esda.getisord import G_Local
from statsmodels.stats.multitest import multipletests
import logging
from datetime import datetime, timezone

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(message)s",
    handlers=[logging.StreamHandler()],
)


def optimize_gistar_bandwidth(
    gdf: gpd.GeoDataFrame,
    variable: str,
    bw_range: np.ndarray = np.arange(1000, 15000, 1000),
    fdr_alpha: float = 0.05,
    min_neighbors: int = 3,
    row_standardize: bool = True,
    seed: int = 42,
) -> dict:
    """Select a Getis-Ord Gi* bandwidth by FDR-stabilized hotspot yield.

    Sweeps `bw_range` (metres), validating weights topology and the
    minimum-neighbor floor at each step, and returns the bandwidth whose
    Gi* output maximizes significant hotspots after FDR correction.
    """
    # 1. CRS enforcement — distance thresholds are only valid in a metre CRS.
    if gdf.crs is None or not gdf.crs.is_projected:
        target_crs = gdf.estimate_utm_crs()
        logging.info("Non-projected CRS; auto-projecting to %s", target_crs)
        gdf = gdf.to_crs(target_crs)
    else:
        target_crs = gdf.crs

    # 2. Deterministic ordering — Gi* is order-sensitive, so sort by a stable id.
    sort_col = "case_id" if "case_id" in gdf.columns else gdf.index.name or None
    if sort_col and sort_col in gdf.columns:
        gdf = gdf.sort_values(sort_col).reset_index(drop=True)

    y = gdf[variable].to_numpy(dtype=float)
    if np.any(~np.isfinite(y)):
        raise ValueError("Analysis variable contains NaN/Inf; impute or filter first.")

    results_log, best = [], None

    # 3. Bandwidth sweep.
    for bw in bw_range:
        try:
            w = libpysal.weights.DistanceBand.from_dataframe(
                gdf, threshold=float(bw), binary=False, silence_warnings=True
            )
            if row_standardize:
                w.transform = "R"

            if w.n != len(gdf):
                logging.warning("BW %dm: topology mismatch (w.n=%d); skipping.", bw, w.n)
                continue

            # cardinalities is a dict keyed by feature id -> neighbor count.
            neighbor_counts = np.array(list(w.cardinalities.values()))
            if neighbor_counts.min() < min_neighbors:
                logging.info("BW %dm: below min_neighbors floor; skipping.", bw)
                continue

            gi = G_Local(y, w, star=True, seed=seed)  # star=True => Gi*

            # FDR — multipletests -> (reject, pvals_corrected, alphacSidak, alphacBonf)
            _, p_fdr, _, _ = multipletests(gi.p_sim, alpha=fdr_alpha, method="fdr_bh")

            sig_mask = p_fdr < fdr_alpha
            hotspot_count = int(sig_mask.sum())

            # Penalize hotspots resting on the neighbor floor (fragmentation proxy).
            fragile = int(np.sum(neighbor_counts[sig_mask] <= min_neighbors))
            score = hotspot_count - 0.1 * fragile

            results_log.append({
                "bandwidth_m": int(bw),
                "hotspot_count": hotspot_count,
                "fragile_hotspots": fragile,
                "score": float(score),
                "min_neighbors_obs": int(neighbor_counts.min()),
                "fdr_alpha": fdr_alpha,
                "timestamp": datetime.now(timezone.utc).isoformat(),
            })

            if best is None or score > best["score"]:
                best = {
                    "score": float(score),
                    "bandwidth_m": int(bw),
                    "weights": w,
                    "result": {
                        "z_scores": gi.Zs,        # G_Local exposes z-scores as .Zs
                        "p_values": gi.p_sim,
                        "p_fdr": p_fdr,
                        "hotspot_flag": sig_mask,
                    },
                }

        except Exception as exc:  # noqa: BLE001 — log and continue the sweep
            logging.error("BW %dm failed: %s", bw, exc)
            continue

    if best is None:
        raise RuntimeError(
            "No valid bandwidth in range. Widen bw_range or lower min_neighbors."
        )

    logging.info("Optimal bandwidth: %dm (score %.2f)", best["bandwidth_m"], best["score"])
    return {
        "optimal_bandwidth_m": best["bandwidth_m"],
        "weights_matrix": best["weights"],
        "results": best["result"],
        "optimization_log": pd.DataFrame(results_log),
        "crs_used": target_crs.to_epsg(),
        "fdr_alpha": fdr_alpha,
        "seed": seed,
    }

The returned optimization_log is the artifact that makes the choice reviewable: one row per candidate bandwidth with its significant-hotspot count, fragile-hotspot penalty, and the observed minimum neighbor count. A reviewer can read the peak straight off that table instead of trusting a single number, and the table doubles as the response curve plotted in the diagram above.

Validation & Edge Cases

Three failure modes account for nearly every wrong bandwidth in practice. Each has a diagnostic you can assert on before publishing the heatmap.

1. The whole range is skipped — a flat, empty log. If every bandwidth trips the min_neighbors floor, optimize_gistar_bandwidth raises and results_log is empty. This is the signature of a sweep whose lower bound is far below the spatial scale of the layer (e.g. 1 km bands on tracts that are 10 km apart). Inspect the nearest-neighbor distribution to set a sane floor before sweeping:

from libpysal.weights import KNN
nn = KNN.from_dataframe(gdf, k=1)
d1 = np.array([d for d in nn.weights.values()]).ravel()  # 1st-NN distances (m)
logging.info("1st-NN distance: median=%.0fm, 95th pct=%.0fm", np.median(d1), np.percentile(d1, 95))
# Start bw_range at roughly the 95th-percentile 1st-NN distance so rural units clear the floor.

2. A flat-topped score with no clear peak. If several adjacent bandwidths return near-identical hotspot_count, the response curve is being read at too coarse a step. Re-sweep at finer resolution around the plateau rather than accepting the first tie:

log = result["optimization_log"]
top = log.nlargest(3, "score")["bandwidth_m"].to_numpy()
if top.max() - top.min() <= log["bandwidth_m"].diff().median():
    fine = np.arange(top.min() - 500, top.max() + 600, 250)
    logging.info("Flat peak near %s; refining sweep over %s", sorted(top), fine.tolist())
    # result = optimize_gistar_bandwidth(gdf, variable, bw_range=fine, ...)

3. Counts masquerading as a hotspot field. If the chosen bandwidth lights up exactly the densest, most populous units, you are almost certainly mapping population, not risk. Check the rank correlation between the analysis variable and population — a near-perfect correlation means the input is an unadjusted count:

from scipy.stats import spearmanr
rho, _ = spearmanr(gdf[variable], gdf["population"])
if abs(rho) > 0.9:
    logging.error("Input tracks population (rho=%.2f); use an age-standardized rate.", rho)

For very large layers (N > 50k units) the dense DistanceBand matrix dominates memory at wide bandwidths. Hold the weights sparse, cap the upper end of bw_range, and where adaptive neighborhoods are acceptable substitute libpysal.weights.KNN with a fixed k derived from the winning fixed distance, or shard by jurisdiction rather than loading a national layer at once.

Compliance Notes

The bandwidth is a parameter that changes which communities get flagged, so its provenance has to be logged. For regulatory defensibility, persist the following to a JSON audit trail alongside the exported heatmap layer:

• The full optimization_log — every candidate bandwidth and its FDR-significant count — so a reviewer can confirm the peak was selected, not hand-picked.
• optimal_bandwidth_m and the weights scheme actually used (distance band vs KNN fallback), since this defines every downstream z-score.
• CRS EPSG code the sweep ran in, plus confirmation the input was reprojected from its source CRS — a bandwidth computed in the wrong units is silently wrong.
• FDR method and α, plus the raw-vs-corrected significant counts at the chosen bandwidth, so the number of suppressed false hotspots is visible.
• The random seed used for Gi* permutation inference — pseudo p-values are stochastic and irreproducible without it.
• Library versions and a SHA-256 hash of the input layer to lock data lineage, and confirmation that the de-identification gate ran before the layer reached this stage.
• Output schema: export a validated GeoDataFrame with explicit columns — gi_z, gi_p_raw, gi_p_fdr, hotspot_flag, bandwidth_m — and attach ISO 19115-compliant metadata before any interagency handoff.

Frequently Asked Questions

Should I optimize a fixed distance band or use adaptive KNN bandwidths? Use a fixed distance band when transmission risk is genuinely distance-driven and unit sizes are comparable, because it preserves true spatial continuity. Switch to adaptive KNN when population density is so heterogeneous that no single distance clears the minimum-neighbor floor everywhere — KNN guarantees neighbor count at the cost of distorting distance, so log which scheme ran.

Why score on FDR-significant hotspot count instead of maximizing the Gi z-scores?* Wider neighborhoods can inflate individual z-scores while erasing real local structure, so raw z is a misleading objective. FDR-significant count rewards bandwidths that find many genuinely defensible hotspots and penalizes both the noisy under-connected regime and the over-smoothed regime — it tracks the single peak of the response curve.

Can I reuse one bandwidth across a multi-week surveillance run? Re-optimize whenever the underlying geography or denominator changes materially; a bandwidth tuned on one week’s incidence layer is not guaranteed optimal on the next. If you must hold it fixed for comparability, log that decision and the date of the last optimization in the audit trail.

What if two bandwidths tie on significant hotspot count? Prefer the smaller bandwidth that has fewer hotspots resting on the minimum-neighbor floor — that is exactly what the fragility penalty in the score encodes. If the tie persists, refine the sweep at a finer step around the plateau as shown in the validation section rather than picking arbitrarily.

Related Topics

• Getis-Ord Gi* Hotspot Detection — the parent method family this bandwidth feeds, with the shared weights and FDR conventions.
• Calculating Local Moran’s I for Infectious Disease Outbreaks — the sibling local statistic when you need HH/LL cluster classification rather than distance-based hotspot intensity.
• Coordinate Reference Systems for Public Health — the metre-based CRS enforcement every distance threshold here depends on.