Choosing Optimal Bin Sizes for Urban Mobility Heatmaps

import numpy as np
import geopandas as gpd
import pandas as pd
from shapely.geometry import box
from scipy.spatial import KDTree


def compute_optimal_bin_width(
    gdf: gpd.GeoDataFrame,
    percentile: float = 0.75,
    min_width: float = 25.0,
    max_width: float = 1000.0,
) -> float:
    """Derive spatial bin width from the nearest-neighbour distance distribution.

    Parameters
    ----------
    gdf:        GeoDataFrame with point geometry; any CRS accepted.
    percentile: Fraction of the NN distribution to use as the bin width (0–1).
    min_width:  Hard floor in metres to prevent degenerate micro-bins.
    max_width:  Hard ceiling in metres to prevent computationally expensive mega-grids.

    Returns
    -------
    Bin width in metres (float), clamped to [min_width, max_width].

    Raises
    ------
    ValueError: if gdf is empty or has fewer than 2 points.
    """
    if gdf is None or len(gdf) < 2:
        raise ValueError("GeoDataFrame must contain at least 2 points.")

    # Always project to metric CRS — never compute distances in EPSG:4326
    if gdf.crs is None or gdf.crs.is_geographic:
        gdf = gdf.to_crs("EPSG:3857")

    coords = np.column_stack((gdf.geometry.x, gdf.geometry.y))
    tree = KDTree(coords)
    # k=2: index 0 is self (distance 0), index 1 is true nearest neighbour
    distances, _ = tree.query(coords, k=2)
    nn_dists = distances[:, 1]

    bin_width = float(np.percentile(nn_dists, percentile * 100))
    return float(np.clip(bin_width, min_width, max_width))


def build_mobility_heatmap(
    gdf: gpd.GeoDataFrame,
    bin_width: float,
    time_col: str = "timestamp",
    window_min: int = 15,
) -> gpd.GeoDataFrame:
    """Construct a regular spatial grid and aggregate point counts per time bin.

    Parameters
    ----------
    gdf:        Projected GeoDataFrame with point geometry and a timestamp column.
    bin_width:  Spatial bin side length in metres (use compute_optimal_bin_width).
    time_col:   Name of the datetime column in gdf.
    window_min: Temporal bin width in minutes (floor-aligned).

    Returns
    -------
    GeoDataFrame with columns [cell_id, time_bin, count, geometry].
    Empty cells (count == 0) are excluded to avoid zero-inflation.

    Raises
    ------
    ValueError: if gdf is empty, bin_width <= 0, or time_col is missing.
    """
    if gdf is None or gdf.empty:
        raise ValueError("GeoDataFrame is empty.")
    if bin_width <= 0:
        raise ValueError(f"bin_width must be positive, got {bin_width}.")
    if time_col not in gdf.columns:
        raise ValueError(f"Column '{time_col}' not found in GeoDataFrame.")

    # Ensure metric projection
    if gdf.crs is None or gdf.crs.is_geographic:
        gdf = gdf.to_crs("EPSG:3857")

    minx, miny, maxx, maxy = gdf.total_bounds
    x_edges = np.arange(minx, maxx + bin_width, bin_width)
    y_edges = np.arange(miny, maxy + bin_width, bin_width)

    # Vectorised grid generation
    cells = [
        box(x, y, x + bin_width, y + bin_width)
        for x in x_edges
        for y in y_edges
    ]
    grid = gpd.GeoDataFrame({"geometry": cells}, crs=gdf.crs)
    grid["cell_id"] = grid.index

    # Spatial join — drop points that fall outside the bounding box extension
    joined = gpd.sjoin(gdf, grid, how="inner", predicate="within")

    # Floor-align timestamps to temporal bin boundaries
    joined["time_bin"] = pd.to_datetime(joined[time_col]).dt.floor(f"{window_min}min")

    # Aggregate and merge geometry back
    counts = (
        joined.groupby(["cell_id", "time_bin"])
        .size()
        .reset_index(name="count")
    )
    heatmap = counts.merge(grid[["cell_id", "geometry"]], on="cell_id")

    # Exclude empty cells — retaining zeros distorts normalisation
    return gpd.GeoDataFrame(
        heatmap[heatmap["count"] > 0].reset_index(drop=True),
        crs=gdf.crs,
    )


# --- Usage ---
# optimal_w = compute_optimal_bin_width(mobility_gdf, percentile=0.75)
# heatmap = build_mobility_heatmap(mobility_gdf, bin_width=optimal_w, window_min=15)

from libpysal.weights import Queen
from esda.moran import Moran

# 1. Shape sanity
assert heatmap.columns.tolist() == ["cell_id", "time_bin", "count", "geometry"]
assert (heatmap["count"] > 0).all(), "Zero-count cells leaked through — filter before analysis."

# 2. Coefficient of Variation — target < 0.8 per time slice
for t_bin, group in heatmap.groupby("time_bin"):
    cv = group["count"].std() / group["count"].mean()
    if cv >= 0.8:
        print(f"WARNING: CV={cv:.2f} at {t_bin} — bins may be too fine (sparse) or too coarse (saturated).")

# 3. Spatial autocorrelation (Moran's I) — target < 0.3 for baseline heatmaps
# Run on the densest time slice to capture worst-case clustering
peak_slice = heatmap.loc[heatmap.groupby("time_bin")["count"].transform("sum").idxmax()]
w = Queen.from_dataframe(gpd.GeoDataFrame(peak_slice, geometry="geometry"))
w.transform = "r"
moran = Moran(peak_slice["count"].values, w)
print(f"Moran's I = {moran.I:.3f}  (p = {moran.p_sim:.3f})")
if moran.I >= 0.3:
    print("Spatial autocorrelation is high — consider halving bin_width and re-validating.")