Selecting KDE Bandwidth for Synthetic Heatmaps

When a kernel density heatmap built from synthetic points collapses two genuinely distinct hotspots into one smooth blob, the culprit is almost never the kernel shape — it is a bandwidth chosen by a rule of thumb that assumes your data is a single Gaussian mound.

Part of Density Mapping & Heat Generation: that workflow turns generated point sets into continuous density surfaces, and the bandwidth is the single parameter that decides whether those surfaces preserve the clustering structure the points were sampled to represent or wash it out. This page treats bandwidth as the bias/variance control it actually is, compares the closed-form Scott and Silverman rules against likelihood cross-validation, and shows why an over-smoothed heatmap silently destroys spatial signal that every downstream consumer then trusts.

Root Cause: Bandwidth Is the Bias/Variance Knob

A two-dimensional kernel density estimate reconstructs a continuous intensity surface from $n$ sample points $\mathbf{x}_i$ by placing a scaled kernel on each and summing:

$$\hat{f}_h(\mathbf{x}) = \frac{1}{n\,h^2} \sum_{i=1}^{n} K\!\left(\frac{\lVert \mathbf{x} - \mathbf{x}_i \rVert}{h}\right)$$

The kernel $K$ (Gaussian, Epanechnikov) barely matters for the result. The bandwidth $h$ — the radius over which each point’s mass is spread — controls everything. It sits directly on the bias/variance trade-off. When $h$ is too small the estimate has low bias but high variance: every sampling accident becomes its own spurious peak, and the heatmap is a spiky field of individual points rather than a density. When $h$ is too large the estimate has low variance but high bias: mass bleeds across real gaps, adjacent clusters merge, and the surface converges toward a single broad hump that no longer encodes where the points actually concentrate.

The over-smoothing failure is the dangerous one because it looks clean. A smooth, plausible heatmap reads as correct, so nobody questions it — yet if the bandwidth exceeds the separation between two clusters, the trough between them fills in and the two modes fuse into one. The heatmap now asserts a spatial structure that contradicts the point set it was built from.

Scott’s and Silverman’s rules are the usual defaults, and both are derived by minimizing the asymptotic mean integrated squared error under the assumption that the underlying density is Gaussian. For a $d$-dimensional sample Scott’s rule sets $h \propto n^{-1/(d+4)}\,\hat{\sigma}$, so in the plane ($d=2$) the bandwidth shrinks only as $n^{-1/6}$ and scales with the marginal standard deviation. That is exactly wrong for the clustered, multimodal point fields that synthetic spatial generation produces: the sample standard deviation is inflated by the distance between clusters, so the rule prescribes a bandwidth wide enough to smear across them. The rule optimizes for a shape your data does not have. Likelihood cross-validation makes no distributional assumption — it selects the $h$ that best predicts held-out points — which is why it recovers separated modes that the closed-form rules erase.

Minimal Reproducer: Watch Scott’s Rule Merge Two Clusters

Generate two well-separated clusters and let SciPy’s default (Scott’s rule) pick the bandwidth. Because the marginal spread is dominated by the gap between clusters, the rule returns a bandwidth wider than the gap itself, and the estimated surface has a single peak where the data has two.

pythonCopy
# requirements.txt — pin majors; let patch releases float
numpy==1.26.*
scipy==1.11.*
scikit-learn==1.4.*
geopandas==0.14.*
pyproj==3.*

pythonCopy
import numpy as np
from scipy.stats import gaussian_kde

rng = np.random.default_rng(20260701)          # explicit integer seed, never wall-clock
# Two clusters ~600 m apart, in a metric CRS (metres), tight 40 m spread.
a = rng.normal([0.0,   0.0], 40.0, size=(400, 2))
b = rng.normal([600.0, 0.0], 40.0, size=(400, 2))
pts = np.vstack([a, b]).T                       # shape (2, N) for scipy

kde = gaussian_kde(pts)                          # Scott's rule by default
h_scott = kde.factor * pts.std(axis=1).mean()    # effective bandwidth in metres
print(f"Scott bandwidth ~ {h_scott:6.1f} m  (cluster gap = 600 m)")

# Sample the density along the x-axis through both cluster centres.
xs = np.linspace(-200, 800, 200)
line = np.vstack([xs, np.zeros_like(xs)])
dens = kde(line)
n_modes = int(((dens[1:-1] > dens[:-2]) & (dens[1:-1] > dens[2:])).sum())
print(f"modes recovered along axis: {n_modes}")   # -> 1, the two clusters have merged

A bandwidth of a few hundred metres against a 600 m gap is enough to fill the trough. The rule is not broken; it is answering the wrong question, because the sample is not one Gaussian.

Fix: Reproject to Metres, Then Cross-Validate the Bandwidth

Two changes recover the structure. First, do the estimate in a metric or equal-area CRS so the bandwidth is an isotropic distance — a degree of longitude is not a degree of latitude away from the equator, and a bandwidth in degrees stretches the kernel into an ellipse. Second, replace the closed-form rule with likelihood cross-validation: search a grid of bandwidths and keep the one that maximizes the leave-one-out log-likelihood of held-out points. This is exactly what GridSearchCV does when it scores KernelDensity by its score (total log-likelihood), and it makes no assumption about the density’s shape.

pythonCopy
import numpy as np
import geopandas as gpd
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV

def project_to_metric(gdf: gpd.GeoDataFrame) -> np.ndarray:
    """Reproject points into an equal-area metric CRS so bandwidth is isotropic metres."""
    if gdf.crs is None:
        raise ValueError("point layer has no CRS; refuse to guess")
    metric = gdf.to_crs("EPSG:6933")             # World Cylindrical Equal Area, metres
    return np.column_stack([metric.geometry.x, metric.geometry.y])

def select_bandwidth(xy: np.ndarray, seed: int = 20260701) -> float:
    """Likelihood cross-validation: maximize held-out log-likelihood over a log grid."""
    # Search around a data-driven scale: median nearest-neighbour distance to typical spread.
    lo, hi = np.log10(5.0), np.log10(400.0)      # metres; bracket the plausible range
    grid = {"bandwidth": np.logspace(lo, hi, 40)}
    search = GridSearchCV(
        KernelDensity(kernel="gaussian"),
        grid,
        cv=5,                                     # 5-fold ~ leave-one-out for held-out LL
        n_jobs=-1,
    )
    search.fit(xy)
    h = float(search.best_params_["bandwidth"])
    if np.isclose(h, grid["bandwidth"][0]) or np.isclose(h, grid["bandwidth"][-1]):
        raise ValueError(f"selected h={h:.1f} m sits on the search boundary; widen the grid")
    return h

def density_grid(xy: np.ndarray, h: float, cell: float = 20.0) -> np.ndarray:
    """Evaluate the CV-selected KDE on a raster grid for the heatmap."""
    kde = KernelDensity(kernel="gaussian", bandwidth=h).fit(xy)
    xmin, ymin = xy.min(axis=0) - 3 * h
    xmax, ymax = xy.max(axis=0) + 3 * h
    gx, gy = np.meshgrid(np.arange(xmin, xmax, cell), np.arange(ymin, ymax, cell))
    flat = np.column_stack([gx.ravel(), gy.ravel()])
    return np.exp(kde.score_samples(flat)).reshape(gx.shape)   # density surface

The logspace grid searches bandwidth on a log scale because the log-likelihood curve is roughly symmetric in $\log h$, not in $h$. The boundary check is not optional: if the optimum lands on the smallest or largest candidate, cross-validation has not actually found a maximum, and you must widen the grid before trusting the result. The scale of the surface here is the same problem that scaling density-based spatial generation with Dask addresses once the grid no longer fits in memory — bandwidth selection runs on a sample, then the chosen $h$ is applied tile-by-tile.

Verification Step: Assert Structure Survives

Gate the bandwidth in CI. The two assertions that matter are that cross-validation found an interior optimum and that the resulting surface still resolves the number of modes the generator put in. A peak-count along the density axis is a cheap, direct proxy for “did over-smoothing erase a mode.”

pythonCopy
import numpy as np

def test_cv_bandwidth_preserves_structure(xy, expected_modes: int):
    h = select_bandwidth(xy)                     # raises if optimum is on the grid edge

    # 1. Cross-validated bandwidth must be well below the inter-cluster separation.
    #    Estimate separation as the span of cluster-centre coordinates on the long axis.
    span = xy[:, 0].ptp()
    assert h < 0.25 * span, f"h={h:.1f} m too wide for span {span:.1f} m -> merging risk"

    # 2. The density surface must recover the expected number of modes.
    surf = density_grid(xy, h)
    row = surf[surf.shape[0] // 2]               # slice through the cluster axis
    peaks = int(((row[1:-1] > row[:-2]) & (row[1:-1] > row[2:])).sum())
    assert peaks >= expected_modes, f"recovered {peaks} modes, expected {expected_modes}"

    # 3. Reproducibility: same seed -> same bandwidth, bitwise.
    assert select_bandwidth(xy) == h, "bandwidth selection is non-deterministic"

Pair the mode-count gate with the quantitative distribution checks in realism metrics evaluation: a heatmap that passes the peak count but shifts the density mass will still fail a Wasserstein comparison against the source intensity, so both gates run together.

Edge Cases & Gotchas

Bandwidth in degrees on an unprojected layer. Fitting a KDE directly on EPSG:4326 coordinates makes the bandwidth a longitude/latitude quantity, and one degree of longitude collapses toward the poles. At 60° latitude a “square” kernel is twice as wide east–west as north–south in metres, so the heatmap smears along parallels. Always reproject to a metric or equal-area CRS (EPSG:6933, or a local UTM zone) before selecting or applying $h$, and carry the CRS as part of the artifact.

Coincident and duplicated points drive the LSCV bandwidth to zero. If the synthetic generator emits exact duplicates — common when points are snapped to a grid or an H3 cell centroid — the leave-one-out likelihood can be maximized by an arbitrarily small bandwidth that places infinite density on the repeats. Cross-validation then selects the smallest candidate and trips the boundary guard. De-duplicate, or jitter coincident points by a sub-metre amount, before bandwidth selection. The same discreteness that helps generating urban point patterns using Poisson processes reproducible can produce these ties.

Edge bias near the study-area boundary. A fixed-bandwidth KDE underestimates density within one bandwidth of the domain edge, because half the kernel falls outside the observed region with no points to fill it. Hotspots that hug a coastline or an administrative border read as artificially cool. Correct with a boundary-reflection or a normalization by the fraction of each kernel inside the domain, and be aware this interacts with the edge treatment discussed for point patterns.

One global bandwidth cannot serve variable density. A single $h$ optimal for a dense downtown over-smooths it while still being too wide for a sparse rural fringe, or vice versa. When the intensity varies by orders of magnitude across the map, a fixed bandwidth is the wrong model — move to an adaptive (variable-bandwidth) KDE that shrinks $h$ in dense regions, and cross-validate the pilot bandwidth rather than a single global value.

Frequently Asked Questions

Related

• Density Mapping & Heat Generation — the parent workflow that turns generated point sets into the continuous surfaces this page tunes.
• Scaling Density-Based Spatial Generation with Dask — apply the chosen bandwidth tile-by-tile once the grid no longer fits in memory.
• Generating Urban Point Patterns Using Poisson Processes — the point sets that feed the heatmap, and the source of coincident-point ties.
• Realism Metrics Evaluation — score the resulting density surface against the source intensity so mode count and mass are both checked.
• Memory-Optimized Large-Grid Spatial Generation — evaluate the tuned kernel over rasters too large to hold in RAM.