Kernel Density Surfaces

A kernel density surface is a continuous estimate of mobility intensity, built by centring a smooth kernel on every observed point or trajectory vertex and summing their overlapping contributions across a regular grid to produce a heat surface free of the blocky artefacts that plague count-based binning.

Prerequisites Checklist

Kernel density estimation sits near the end of a spatial aggregation workflow: it consumes cleaned, projected coordinates and produces a raster that other tools consume. Before building a surface, confirm the following. This topic is one branch of Spatial Indexing & Density Aggregation, which frames how continuous surfaces, discrete grids, and flow matrices complement one another.

• Python ≥ 3.10 with numpy ≥ 1.24 and either scipy ≥ 1.11 (for scipy.stats.gaussian_kde) or scikit-learn ≥ 1.3 (for sklearn.neighbors.KernelDensity).
• pyproj ≥ 3.6 for datum transformations, imported with always_xy=True to avoid axis-order bugs.
• rasterio ≥ 1.3 to attach an affine transform and CRS to the output array, and shapely ≥ 2.0 if you derive contour polygons.
• Clean input. Points should already pass through GPS filtering — raw multipath jumps inflate density in the wrong places. See GPS precision and error handling for the upstream cleaning contract.
• Input schema: a table with lon (float64), lat (float64) in EPSG:4326, optionally weight (float64) for intensity-weighted surfaces (dwell seconds, passenger counts), and for time-sliced work a ts (UTC datetime64) column.
• Projection context: GPS emits angular EPSG:4326. Bandwidth, cell size, and every distance in this pipeline must be expressed in a metric projected CRS — this is flagged explicitly in every code block below.

Error & Problem Taxonomy

Most failed density surfaces are not bugs in the estimator; they are misconfigurations of scale, units, or geometry. Classifying the symptom tells you which knob to turn — an oversmoothed blob and an aliased staircase have opposite fixes, and confusing them wastes tuning cycles.

Deterministic Pipeline Overview

The five-stage sequence below is the standard production order. Each stage is fully vectorized — no Python-level row iteration — and each records the parameters it used so a surface can be reproduced and audited later.

1. Project to a metric CRS — reproject lon/lat into a UTM zone (or a national grid) derived from the data centroid, so every subsequent length is in metres.
2. Choose kernel and bandwidth — pick a Gaussian kernel and a bandwidth h, either from a rule-of-thumb estimator (Scott, Silverman) as a starting point or set explicitly from the phenomenon scale.
3. Evaluate on a raster grid — construct a regular grid whose extent is the buffered data envelope and whose cell size is a fraction of h; evaluate the density estimator at every cell centre.
4. Normalise the surface — scale the array to integrate to one over the area (a probability density), or convert to interpretable units such as points per square kilometre.
5. Export raster and contours — pair the array with an affine transform and CRS and write a GeoTIFF; optionally trace iso-density contour polygons for vector overlays and cartographic labelling.

The concrete, runnable form of these five stages — including the rasterio affine transform and the singular-covariance edge cases — is worked through end to end in generating kernel density heatmaps from mobility data.

Implementation Walkthrough

The function below fits a Gaussian KDE on projected metric coordinates and evaluates it over a raster grid, returning both the density array and the affine transform needed to georeference it. Distances, bandwidth, and cell size are all in metres — never raw EPSG:4326 degrees. Adjust the UTM zone to your data’s extent, or pass an explicit target_crs_epsg.

from __future__ import annotations

import logging
from dataclasses import dataclass

import numpy as np
import pandas as pd
from pyproj import Transformer
from scipy.stats import gaussian_kde

logger = logging.getLogger(__name__)


@dataclass
class DensitySurface:
    """Georeferenced kernel density surface."""
    density: np.ndarray          # 2D array, north-up (row 0 = top / max y)
    transform: tuple[float, ...]  # affine (a, b, c, d, e, f) in the metric CRS
    crs_epsg: int                 # EPSG code of the metric CRS the grid lives in
    bandwidth_m: float            # effective bandwidth in metres
    cell_size_m: float            # grid resolution in metres


def kernel_density_surface(
    df: pd.DataFrame,
    lon_col: str = "lon",
    lat_col: str = "lat",
    weight_col: str | None = None,
    target_crs_epsg: int | None = None,
    bandwidth_m: float | None = None,
    cell_size_m: float = 50.0,
    pad_bandwidths: float = 3.0,
) -> DensitySurface:
    """
    Estimate a continuous kernel density surface from mobility points.

    Parameters
    ----------
    df : DataFrame with lon_col, lat_col in EPSG:4326 and optional weight_col.
    weight_col : per-point intensity (dwell seconds, ridership). None = uniform.
    target_crs_epsg : metric CRS for evaluation. NEVER 4326 — bandwidth would
                      be interpreted in degrees. If None, a UTM zone is derived
                      from the data centroid.
    bandwidth_m : Gaussian bandwidth in METRES. If None, Scott's rule is used
                  (data-driven, but tends to oversmooth clustered mobility data).
    cell_size_m : raster cell size in metres; keep it below the bandwidth.
    pad_bandwidths : buffer the grid envelope by this many bandwidths so the
                     kernel tails are not clipped at the study-area edge.

    Returns
    -------
    DensitySurface with a north-up density array and an affine transform.

    Raises
    ------
    ValueError : missing columns, empty frame, or fewer than 3 usable points.
    """
    required = {lon_col, lat_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Input DataFrame missing columns: {missing}")

    pts = df[[lon_col, lat_col]].dropna()
    if pts.empty:
        raise ValueError("No valid coordinates after dropping NaN rows.")
    if len(pts) < 3:
        # gaussian_kde needs >= 2 points and a non-singular covariance;
        # 3 distinct points is the practical floor for a 2D estimate.
        raise ValueError(f"Need >= 3 points for a 2D KDE, got {len(pts)}.")

    # ── Stage 1: project to a metric CRS ─────────────────────────────────────
    # Derive a UTM zone from the centroid unless the caller pinned a CRS.
    # NEVER fit the KDE on EPSG:4326 — a degree is not a metre and the
    # bandwidth would smear anisotropically with latitude.
    lon_c, lat_c = float(pts[lon_col].mean()), float(pts[lat_col].mean())
    if target_crs_epsg is None:
        zone = int((lon_c + 180) // 6) + 1
        target_crs_epsg = (32600 if lat_c >= 0 else 32700) + zone

    to_metric = Transformer.from_crs("EPSG:4326", f"EPSG:{target_crs_epsg}", always_xy=True)
    x, y = to_metric.transform(pts[lon_col].to_numpy(), pts[lat_col].to_numpy())
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    # Guard against a singular covariance (all points identical / colinear).
    if np.ptp(x) < 1e-6 and np.ptp(y) < 1e-6:
        raise ValueError("All points are effectively identical; covariance is singular.")

    weights = None
    if weight_col is not None:
        weights = df.loc[pts.index, weight_col].to_numpy(dtype=float)
        weights = np.where(np.isfinite(weights) & (weights > 0), weights, 0.0)
        if weights.sum() == 0:
            raise ValueError("All weights are zero or invalid.")

    # ── Stage 2: kernel + bandwidth ──────────────────────────────────────────
    coords = np.vstack([x, y])
    kde = gaussian_kde(coords, weights=weights)  # Scott's rule by default
    if bandwidth_m is not None:
        # gaussian_kde scales its factor by the data covariance, so convert a
        # metric bandwidth into the estimator's dimensionless factor.
        data_std = np.sqrt(np.mean([x.var(ddof=1), y.var(ddof=1)]))
        kde.set_bandwidth(bandwidth_m / data_std)
    effective_bw = float(np.sqrt(kde.covariance[0, 0]))

    # ── Stage 3: evaluate on a raster grid ───────────────────────────────────
    pad = pad_bandwidths * effective_bw
    x_min, x_max = x.min() - pad, x.max() + pad
    y_min, y_max = y.min() - pad, y.max() + pad
    n_cols = max(1, int(np.ceil((x_max - x_min) / cell_size_m)))
    n_rows = max(1, int(np.ceil((y_max - y_min) / cell_size_m)))

    xs = x_min + (np.arange(n_cols) + 0.5) * cell_size_m
    ys = y_min + (np.arange(n_rows) + 0.5) * cell_size_m
    gx, gy = np.meshgrid(xs, ys)
    grid = np.vstack([gx.ravel(), gy.ravel()])
    density = kde(grid).reshape(n_rows, n_cols)

    # North-up: row 0 should map to the TOP (max y), so flip vertically.
    density = np.flipud(density)

    # ── Stage 4: normalise so the surface integrates to ~1 over the area ──────
    cell_area = cell_size_m * cell_size_m
    total = density.sum() * cell_area
    if total > 0:
        density = density / total

    # Affine transform (GDAL order): pixel (col,row) -> metric (x,y), north-up.
    transform = (cell_size_m, 0.0, x_min, 0.0, -cell_size_m, y_max)

    logger.info(
        "KDE surface: %d points -> %dx%d grid, bandwidth=%.1f m, cell=%.1f m, CRS=EPSG:%d",
        len(pts), n_rows, n_cols, effective_bw, cell_size_m, target_crs_epsg,
    )
    return DensitySurface(
        density=density,
        transform=transform,
        crs_epsg=target_crs_epsg,
        bandwidth_m=effective_bw,
        cell_size_m=cell_size_m,
    )

The sklearn.neighbors.KernelDensity estimator is the alternative when you need a non-Gaussian kernel (Epanechnikov, tophat) or a tree-backed evaluation that scales to millions of points; it takes an explicit bandwidth argument already expressed in the units of the fitted coordinates, so fitting it on projected metres gives a bandwidth in metres directly, with no covariance conversion.

Geometric & Mathematical Grounding

The estimator. For n observed points at metric locations xᵢ, the kernel density estimate at a location x is f̂(x) = (1 / n·h²) · Σ K((x − xᵢ) / h), where K is the kernel function and h is the bandwidth. Each point contributes a small bump of total mass 1/n; the surface is the sum of all bumps. Because the contributions add linearly, the result is smooth and — unlike a histogram or grid count — independent of where you happen to place the grid origin.

The Gaussian kernel. The most common choice is the isotropic 2D Gaussian, K(u) = (1 / 2π) · exp(−‖u‖² / 2). It has infinite support, so every point technically influences every cell; in practice contributions beyond about three bandwidths are negligible, which is why the pipeline pads the grid envelope by three bandwidths rather than an arbitrary margin. The Epanechnikov kernel is theoretically optimal in a mean-integrated-squared-error sense and has finite support (it is exactly zero beyond h), which makes it cheaper for tree-based evaluation, but the visual difference on mobility heat surfaces is minor.

Bandwidth selection. Scott’s rule sets h = σ · n^(−1/(d+4)) and Silverman’s rule is a close variant; for two dimensions (d = 2) both shrink the bandwidth slowly as n grows. They are derived under the assumption that the underlying density is itself Gaussian, which mobility data rarely is — trips cluster tightly along roads and around venues. As a result these rules systematically oversmooth clustered movement data, so treat the rule-of-thumb value as a floor and tune upward or downward against a held-out validation sample rather than accepting it blindly.

Why the bandwidth must be metric. The bandwidth is a physical smoothing radius that enters the estimator as (x − xᵢ) / h. If x and xᵢ are in degrees, then h is in degrees, and a fixed h corresponds to a different ground distance at every latitude because a degree of longitude contracts toward the poles. Projecting to a metric CRS first — the same discipline that governs every distance and speed computation across this site — makes h = 150 mean 150 metres uniformly, so the kernel is isotropic on the ground and density values are comparable across the whole surface.

Calibration & Parameter Tuning

Bandwidth is the parameter that determines whether a surface is useful or misleading, and the right value depends on both the phenomenon you are mapping and the transport mode that generated the points. The table gives practical starting values; the cell size follows the bandwidth, not the other way around.

Calibrate against your own data rather than these defaults alone. A reliable procedure is to hold out a random 20% of points, build the surface from the remaining 80%, and score the held-out points by the log-likelihood the surface assigns them; sweep the bandwidth and keep the value that maximises held-out likelihood, then nudge it for cartographic taste. For multi-modal datasets, build a separate surface per mode inferred from speed and acceleration profiling rather than one blended surface, because a single bandwidth cannot serve both 40-metre pedestrian dwell clusters and 500-metre highway corridors.

Integration & Compatibility

Kernel density surfaces are one of several aggregation strategies, and they interlock with the others rather than replacing them.

• Discrete Global Grid Systems: grid binning and KDE answer different questions. Discrete grids give exact, joinable counts per cell and are the right tool when you need to attach density to administrative or hexagonal units; KDE gives a smooth, origin-independent surface for cartography and gradient analysis. A common pattern is to compute exact counts on an H3 hexagon grid for tabular reporting and a KDE surface for the accompanying heat map.
• Origin-Destination Flow Matrices: KDE describes where activity concentrates; OD matrices describe how activity moves between zones. Density surfaces of trip endpoints often reveal the natural zone boundaries that make an OD matrix legible.
• Time-sliced surfaces: to animate how a hotspot pulses over a day, build one surface per time bin using a shared grid and bandwidth, then stack them. Choosing the bin width is itself a design decision covered in dynamic time-binning strategies — too fine and each frame is noise, too coarse and the temporal signal is averaged away.
• Coordinate Reference System Mapping: the whole surface is only as trustworthy as its projection. Pin one metric CRS for the entire aggregation stage so surfaces, grids, and contours overlay without re-projection drift.
• Generating kernel density heatmaps from mobility data: the step-by-step deep-dive with the full rasterio export, PNG rendering, and edge-case handling.

FAQ

Why must the bandwidth be expressed in metres rather than degrees?

A kernel bandwidth is a physical smoothing radius. In EPSG:4326 one degree of longitude spans roughly 111 km at the equator but shrinks toward the poles, so a bandwidth of 0.01 degrees means very different ground distances at different latitudes and stretches the kernel anisotropically. Project to a metric CRS such as a UTM zone first, and a bandwidth of 150 means 150 metres everywhere on the grid.

How do I pick a bandwidth for a kernel density surface?

Start from the scale of the phenomenon you want to see. Dwell hotspots around venues and stations resolve well at 50–150 m; corridor-scale flow needs 200–500 m. Rule-of-thumb estimators such as Scott and Silverman give a data-driven starting bandwidth, but they assume a roughly Gaussian point cloud and oversmooth clustered mobility data, so treat them as a lower bound and tune upward against a held-out sample.

What raster cell size should I pair with a given bandwidth?

Set the cell size to roughly one quarter to one half of the bandwidth. Cells larger than the bandwidth alias the smooth surface into blocky steps; cells far smaller than the bandwidth multiply compute and storage cost without adding resolvable detail. A 150 m bandwidth pairs naturally with a 50 m cell.

How is a kernel density surface different from a hexagon or grid count?

Grid binning assigns each point to exactly one discrete cell and counts membership, so the output is a step function whose appearance depends on the arbitrary placement of cell edges. KDE spreads each point over a smooth kernel and sums the overlaps, producing a continuous surface insensitive to grid origin. Use discrete grids for exact aggregation and joins; use KDE for smooth hotspot cartography.

How do I handle edge effects at the boundary of my study area?

A standard kernel places mass outside the study boundary where no data can exist, so density is underestimated within one bandwidth of the edge. Buffer the study area by at least one bandwidth before evaluation and clip afterward, apply an edge-correction weight equal to the reciprocal of the kernel mass falling inside the boundary, or reflect points across the boundary. Document the choice, because it changes absolute density values near coastlines and cordon lines.

Related

• Generating Kernel Density Heatmaps from Mobility Data
• Discrete Global Grid Systems
• Origin-Destination Flow Matrices
• Dynamic Time-Binning Strategies
• GPS Precision & Error Handling

Back to Spatial Indexing & Density Aggregation