Generating kernel density heatmaps from mobility data

Turning a pile of GPS points into a clean, georeferenced heatmap takes four steps: project to a metric CRS → fit a kernel density estimator → evaluate on a raster grid → write a GeoTIFF (and optionally a PNG). The two failures that ruin most first attempts are fitting the estimator on raw EPSG:4326 degrees — which smears the kernel anisotropically with latitude — and forgetting to flip the array so north ends up at the top of the raster. This guide gives a complete, runnable implementation that handles both, plus the empty-frame and singular-covariance edge cases that break naive code in production.

Why this happens

A kernel density estimate places a smooth bump on every point and sums the overlaps into a continuous surface. That surface only means anything if the bump has a consistent physical size, and that requires a metric coordinate system: a bandwidth of 150 has to mean 150 metres at every point on the map, not 0.0013 degrees that shrink toward the poles. The conceptual background — kernel choice, bandwidth theory, edge effects, and the calibration table — lives in the parent guide on kernel density surfaces, itself part of the broader Spatial Indexing & Density Aggregation section.

The other recurring problem is dirty input. A KDE faithfully renders whatever points it is given, so multipath jumps and drift produce density in places nothing actually happened. Run points through the cleaning contract in GPS precision and error handling first; a heatmap of raw pings is a heatmap of receiver noise as much as of movement.

Core pipeline

The four steps below run in strict order — you cannot fit a metric-bandwidth kernel before projecting, and you cannot georeference the raster before you know the grid extent.

Kernel Density Heatmap Pipeline Four sequential stages: Project to Metric CRS, Fit KDE, Evaluate on Grid, and Write GeoTIFF and PNG, connected by arrows. 1. Project to Metric CRS 2. Fit KDE Gaussian, h in m 3. Evaluate on Raster Grid 4. Write GeoTIFF / PNG
  1. Project to a metric CRS. Derive a UTM zone from the data centroid and reproject with pyproj, so bandwidth and cell size are in metres. This is the same projection discipline used across CRS transformation best practices.
  2. Fit the KDE. Fit scipy.stats.gaussian_kde on the projected coordinates, optionally weighted by dwell time or ridership, with the bandwidth set from the phenomenon scale in metres.
  3. Evaluate on a raster grid. Build a regular grid over the buffered data envelope, with cell size below the bandwidth, and evaluate the estimator at every cell centre.
  4. Write the outputs. Attach an affine transform and CRS and write a GeoTIFF with rasterio; optionally normalise to 0–255 and write a PNG for quick sharing.

Production-ready Python implementation

The function below runs the whole pipeline and returns the path to the written GeoTIFF. It projects to a metric CRS, guards the empty-frame, fewer-than-three-points, and all-identical-points cases, evaluates the KDE, flips the array north-up, and writes a georeferenced raster. The pyproj transformer uses always_xy=True so longitude comes before latitude and axis order can never silently swap.

PYTHON 
from __future__ import annotations

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


def mobility_heatmap(
    df: pd.DataFrame,
    out_tif: str,
    lon_col: str = "lon",
    lat_col: str = "lat",
    weight_col: str | None = None,
    bandwidth_m: float = 150.0,
    cell_size_m: float = 50.0,
    pad_bandwidths: float = 3.0,
) -> str:
    """
    Build a georeferenced kernel density heatmap GeoTIFF from GPS points.

    Parameters
    ----------
    df : DataFrame with lon_col, lat_col in EPSG:4326; optional weight_col.
    out_tif : output path for the GeoTIFF.
    weight_col : per-point intensity (dwell seconds, boardings). None = uniform.
    bandwidth_m : Gaussian bandwidth in METRES (physical smoothing radius).
    cell_size_m : raster resolution in metres; keep below the bandwidth.
    pad_bandwidths : pad the grid envelope by this many bandwidths so kernel
                     tails are not clipped at the study-area edge.

    Returns
    -------
    The out_tif path.

    Raises
    ------
    ValueError : missing columns, empty frame, fewer than 3 usable points,
                 or a singular covariance (all points effectively identical).
    """
    required = {lon_col, lat_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")

    pts = df[[lon_col, lat_col]].dropna()
    if pts.empty:
        raise ValueError("Empty input frame: no valid coordinates to map.")
    if len(pts) < 3:
        raise ValueError(f"KDE needs >= 3 points, got {len(pts)}.")

    # ── Step 1: project to a metric CRS (UTM from centroid) ───────────────────
    # NEVER fit the KDE on EPSG:4326 degrees — bandwidth must be metric so the
    # kernel is isotropic on the ground rather than stretched by latitude.
    lon_c, lat_c = float(pts[lon_col].mean()), float(pts[lat_col].mean())
    zone = int((lon_c + 180) // 6) + 1
    metric_epsg = (32600 if lat_c >= 0 else 32700) + zone
    to_metric = Transformer.from_crs("EPSG:4326", f"EPSG:{metric_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)

    # Singular-covariance guard: 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:
        w = df.loc[pts.index, weight_col].to_numpy(dtype=float)
        w = np.where(np.isfinite(w) & (w > 0), w, 0.0)
        if w.sum() == 0:
            raise ValueError("All weights are zero or invalid.")
        weights = w

    # ── Step 2: fit the KDE with a METRIC bandwidth ───────────────────────────
    coords = np.vstack([x, y])
    kde = gaussian_kde(coords, weights=weights)
    data_std = np.sqrt(np.mean([x.var(ddof=1), y.var(ddof=1)]))
    kde.set_bandwidth(bandwidth_m / data_std)  # convert metres -> kde factor

    # ── Step 3: evaluate on a raster grid (cell size < bandwidth) ─────────────
    pad = pad_bandwidths * bandwidth_m
    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)
    density = kde(np.vstack([gx.ravel(), gy.ravel()])).reshape(n_rows, n_cols)

    # North-up: NumPy builds rows from min y upward, but raster row 0 is the
    # TOP (max y). Flip vertically and pair with a negative y-pixel size.
    density = np.flipud(density).astype("float32")

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

    # ── Step 4: write a georeferenced GeoTIFF ─────────────────────────────────
    transform = Affine(cell_size_m, 0.0, x_min, 0.0, -cell_size_m, y_max)
    with rasterio.open(
        out_tif, "w", driver="GTiff",
        height=n_rows, width=n_cols, count=1, dtype="float32",
        crs=f"EPSG:{metric_epsg}", transform=transform, nodata=0.0,
        compress="deflate",
    ) as dst:
        dst.write(density, 1)
        dst.update_tags(bandwidth_m=str(bandwidth_m), cell_size_m=str(cell_size_m))

    return out_tif


def heatmap_to_png(density: np.ndarray, out_png: str) -> str:
    """
    Write a normalised 8-bit PNG for quick visual sharing (NOT georeferenced).
    Scales the density array to 0-255 using its own max as the ceiling.
    """
    import imageio.v3 as iio

    finite_max = np.nanmax(density) if np.isfinite(density).any() else 0.0
    if finite_max <= 0:
        raise ValueError("Density surface is empty or all-zero; nothing to render.")
    scaled = np.clip(density / finite_max, 0.0, 1.0)
    iio.imwrite(out_png, (scaled * 255).astype("uint8"))
    return out_png

Validation block

Run these checks immediately after writing the GeoTIFF. Each targets a distinct failure: an empty or NaN surface, a broken normalisation, or a georeferencing error that would misplace the heatmap on a map.

PYTHON 
import numpy as np
import rasterio


def validate_heatmap(tif_path: str) -> dict:
    """
    Sanity-check a written kernel density GeoTIFF.
    Raises AssertionError on failure; returns a metrics dict on success.
    """
    with rasterio.open(tif_path) as src:
        arr = src.read(1)
        transform = src.transform
        crs = src.crs

    # 1. No NaN or infinities — these break downstream raster algebra.
    assert np.isfinite(arr).all(), "Density surface contains NaN or inf."

    # 2. Density must be non-negative everywhere.
    assert (arr >= 0).all(), "Negative density values found."

    # 3. The surface must integrate to ~1 over the study area.
    #    cell area = |a * e| from the affine transform (metres squared).
    cell_area = abs(transform.a * transform.e)
    integral = float(arr.sum() * cell_area)
    assert 0.95 <= integral <= 1.05, (
        f"Density integral {integral:.3f} not ~1 — check normalisation."
    )

    # 4. Georeferencing: a metric CRS (not EPSG:4326) and a north-up transform.
    assert crs is not None and crs.to_epsg() != 4326, (
        "Raster must be in a metric CRS, not geographic EPSG:4326."
    )
    assert transform.e < 0, (
        "y-pixel size must be negative (north-up); heatmap is likely flipped."
    )

    return {
        "crs_epsg": crs.to_epsg(),
        "shape": arr.shape,
        "density_integral": round(integral, 4),
        "max_density": float(arr.max()),
        "cell_area_m2": round(cell_area, 1),
    }

Key assertions:

  • No NaN / non-negative: a single NaN (often from a bad weight) propagates through raster algebra and corrupts every overlay built on the surface.
  • Integral ≈ 1: confirms the normalisation ran; a value far from one usually means the cell area was omitted from the sum.
  • Metric CRS and negative y-pixel size: catches the two most common georeferencing bugs — leaving the raster in EPSG:4326, and forgetting the vertical flip that puts north at the top.

Common mistakes and gotchas

  • Fitting the KDE on raw longitude/latitude. A bandwidth of 150 in degree-space is meaningless and the kernel stretches anisotropically with latitude. Always project to a metric CRS before fitting — the same rule that governs downstream aggregating GPS pings into H3 hexagon grids and every distance computation.
  • Forgetting the vertical flip. NumPy meshgrids number rows from the minimum y upward, but raster row zero is the north edge. Skip np.flipud (and the negative e term in the affine transform) and your heatmap is mirrored north-to-south even though the coordinates look right.
  • Cell size larger than the bandwidth. This aliases the smooth surface into a blocky staircase. Keep the cell size to roughly a quarter to a half of the bandwidth; a 150 m bandwidth wants a 50 m cell, not a 250 m one.
  • Not padding the grid envelope. If the grid stops exactly at the data extent, the kernel tails are clipped and density is underestimated near the edge. Pad by about three bandwidths before evaluating, then clip to your study area afterward if needed.
  • Ignoring the singular-covariance case. Two points, or many identical points, give gaussian_kde a singular covariance and a LinAlgError. Guard for fewer than three usable points and for a near-zero coordinate spread before fitting, as the implementation above does.
  • Treating the PNG as data. A normalised PNG discards the CRS and affine transform. Share it for a quick look, but keep the GeoTIFF as the analysis-ready artefact for any join, overlay, or measurement.

Back to Kernel Density Surfaces