Spatial Indexing & Density Aggregation

Once a fleet of trajectories has been cleaned, resampled, and segmented, the raw points are still just a scatter of positions—useful for individual inspection but not for the dashboards, exposure models, and traffic estimates that consume movement data. This section covers the aggregation layer: how to convert millions of individual pings into stable, queryable spatial signal without introducing the binning artefacts that quietly corrupt every downstream statistic.

Prerequisites & Scope

This guide assumes Python 3.10+ and comfort with pandas and geopandas, plus the upstream stages that produce the input to everything here. Aggregation only makes sense on data that has already been through ingestion and cleaning in Spatiotemporal Data Foundations & Structures and pattern extraction in Movement Pattern Extraction & Trajectory Analysis. If your points still carry GPS drift or asynchronous clocks, aggregate signal will inherit that error at full strength.

The primary libraries referenced are:

• h3-py ≥ 4.0 — hierarchical hexagonal indexing (H3 v4 API)
• s2sphere — S2 spherical quad-cell indexing
• python-geohash — geohash string encoding for legacy interoperability
• osmnx ≥ 1.9 — downloading and building routable road-network graphs from OpenStreetMap
• leuvenmapmatching / hmmlearn — hidden-Markov-model map matching of noisy traces to edges
• scipy.stats.gaussian_kde and sklearn.neighbors.KernelDensity — continuous density estimation
• geopandas ≥ 0.14 — spatial DataFrames, CRS management, spatial joins
• duckdb ≥ 0.10 with the spatial extension — in-process aggregation over GeoParquet at scale

Every distance, bandwidth, and area computation in this section is performed in a metric projected CRS (a UTM zone or a regional equal-area projection), never raw EPSG:4326. This is stated again wherever coordinates are touched, because degree-space aggregation is the single most common way these pipelines produce plausible-looking but wrong output.

Core Conceptual Model

Aggregation begins with one artefact: a cleaned point stream in which each row is an (entity_id, x, y, t) tuple in a metric CRS, with stay points and moving segments already labelled upstream. From that single stream, three distinct representations branch out, and each produces a different aggregate output.

The grid path: discrete cell indexing

The first branch discretizes continuous coordinates into a discrete global grid system. Each point is assigned a cell key—an H3 hexagon, an S2 quad, or a geohash string—and rows are grouped by that key. The result is a count (or dwell-time sum, or distinct-entity count) per cell: a density map with O(1) joins by cell identifier and a stable hierarchy for zooming between resolutions. Grids ignore topology entirely, which is exactly why they are cheap and why they blur road structure.

The network path: map-matched references

The second branch snaps each noisy point onto the most probable sequence of road edges by map-matching to a road network. Because raw GPS wanders across lane and building boundaries, matching is a probabilistic inference problem, usually solved with a hidden Markov model that balances observation likelihood against transition plausibility. The output is an edge-referenced trace from which you can build link-level traffic volumes, corridor travel times, and—by pairing trip endpoints—an origin-destination flow matrix.

The surface path: continuous density

The third branch estimates a smooth intensity field with kernel density surfaces. Rather than committing to a discrete cell boundary, kernel density estimation (KDE) places a smoothing kernel over every point and sums the contributions onto a metric raster. The result is a continuous surface that avoids the hard edges of a grid but introduces its own critical parameter—the bandwidth, which must be expressed in metres, not degrees.

Architecture Decision Map

The choice of representation, cell system, and output format determines both what questions you can answer and how much compute you spend. The table below maps the decisions you will make first.

Two rules of thumb collapse most of this table. First, if the answer depends on which road was used, you need the network path; otherwise a grid or surface is cheaper and more robust. Second, hexagons are the safer default for movement because their uniform neighbour distance removes the axis bias that squares and quads carry.

Pipeline Integration

This stage sits squarely in the middle of the movement analytics stack. It consumes cleaned, structured trajectories and emits aggregates that dashboards and models treat as inputs.

Upstream — cleaning and structure. Aggregation assumes the input is already correct. Coordinates must be reprojected to a metric CRS, noise filtered, and packaged into consistent trajectory objects. That work happens in Spatiotemporal Data Foundations & Structures. Two dependencies matter especially here. First, unresolved noise from GPS Precision & Error Handling propagates directly into aggregates: a single drifting fix inflates a neighbouring cell’s count and pulls a KDE surface toward an empty street. Second, the schema you inherit from Trajectory Object Design Patterns determines how cleanly you can group by entity and segment—an aggregation is only as reliable as the entity_id and trajectory_id keys it groups on.

Upstream — pattern extraction. Origin-destination matrices depend on trip endpoints, and trip endpoints come from segmentation and stay-point detection in Movement Pattern Extraction & Trajectory Analysis. If moving segments and stops are not separated before aggregation, dwell time at a depot is miscounted as through-traffic and OD flows collapse.

Downstream — dashboards and models. The aggregates produced here are terminal artefacts for many consumers. A binned density map becomes a heatmap tile layer; an OD matrix becomes the input to a gravity or demand model; a KDE raster becomes an exposure covariate. Because these outputs are trusted rather than re-inspected, the validation checks later in this guide are not optional—they are the last gate before wrong numbers reach a decision.

Engineering Pitfalls and Production Gotchas

1. The modifiable areal unit problem (MAUP)

Every discrete grid imposes an arbitrary boundary, and aggregated statistics are sensitive to that choice. The same pings binned at H3 resolution 7 versus 9, or with the tessellation shifted by half a cell, can reorder your hotspots and change the correlation between density and any covariate. The modifiable areal unit problem has two components—a scale effect (cell size) and a zonation effect (cell placement)—and there is no single “correct” resolution. The defensible practice is to compute the aggregation at several resolutions, report how conclusions move, and prefer findings that are stable across them rather than presenting a single grid as ground truth.

2. Grid edge effects

Kernel density estimates and neighbour statistics are biased near the boundary of the study area because cells outside the frame contribute no mass. A hotspot detector will systematically under-report intensity along the edges of your bounding box, sometimes badly enough to hide a genuine peak near the border. Mitigate by buffering the study area, computing density on the buffered extent, and clipping the result back to the region of interest so every reported cell has a full neighbourhood.

3. Resolution mismatch across joined layers

Aggregates are frequently joined—density per cell against demographics per cell, or OD flows against zone attributes. If the two layers use different cell sizes or different tessellations, the join silently redistributes values through an areal-interpolation assumption that is rarely stated. Fix the target resolution once, up front, and reproject or re-bin every layer to it before joining. Never join an H3 resolution-8 layer to a resolution-9 layer without an explicit aggregation step.

4. Bandwidth (or cell size) computed in degrees instead of a metric CRS

This is the highest-frequency correctness bug in density work. A KDE bandwidth of 0.001 or a “cell size” specified as a fraction of a degree is meaningless as a physical distance: one degree of longitude is about 111 km at the equator and shrinks toward the poles, so a degree-space kernel is both wrong in scale and stretched anisotropically with latitude. Always reproject to a metric CRS (a UTM zone or a regional equal-area projection) before choosing a bandwidth, then express it in metres. A 150 m smoothing radius is bw=150, not 0.00135.

5. Hex-versus-square binning bias

Square and quad cells have two distinct neighbour distances—edge-adjacent versus corner-adjacent—so any statistic that walks the neighbourhood (flow, contiguity, diffusion) is biased along the grid axes. Squares also approximate circular kernels poorly, introducing orientation artefacts into density surfaces. Hexagons have a single neighbour distance in all six directions and tile the plane more isotropically, which is why H3 is the usual recommendation for movement aggregation. Choose squares or quads only when you specifically need clean hierarchical nesting, and document that you accepted the axis bias.

Python Tooling Landscape

Validation and Testing Patterns

Aggregation failures are almost always silent: the pipeline runs, produces a plausible map, and quietly drops or double-counts points. Build these checks into your CI/CD data-quality gates so a regression fails the build instead of shipping a wrong dashboard.

Mass conservation. The sum of all per-cell counts must equal the number of input points that fall inside the study region. If it does not, points were dropped—typically by a cell assignment that returned null for out-of-bounds coordinates, or by a join that lost rows.

import geopandas as gpd
import h3
import pandas as pd


def bin_to_h3(gdf: gpd.GeoDataFrame, resolution: int) -> pd.DataFrame:
    """Assign each point an H3 cell and return per-cell counts.

    Coordinates are read as lat/lon from a WGS84 copy for cell assignment,
    but any distance work elsewhere in the pipeline uses a metric CRS.

    Raises:
        ValueError: if the frame is empty or lacks a geometry column.
    """
    if gdf.empty:
        raise ValueError("Empty GeoDataFrame: nothing to bin.")
    if gdf.geometry.name not in gdf.columns:
        raise ValueError("GeoDataFrame has no active geometry column.")

    # H3 indexing operates on geographic lat/lon; take an explicit WGS84 copy
    # so the caller's working CRS (a metric UTM zone) is never mutated.
    wgs = gdf.to_crs(4326)
    cells = [
        h3.latlng_to_cell(pt.y, pt.x, resolution)
        for pt in wgs.geometry
    ]
    counts = pd.Series(cells, name="cell").value_counts().rename_axis("cell")
    return counts.reset_index(name="count")


def assert_mass_conservation(gdf: gpd.GeoDataFrame, binned: pd.DataFrame) -> None:
    """Assert that binned counts sum to the input point count."""
    total_in = len(gdf)
    total_binned = int(binned["count"].sum())
    if total_in != total_binned:
        raise AssertionError(
            f"Mass not conserved: {total_in} input points but "
            f"{total_binned} binned ({total_in - total_binned} lost)."
        )

H3 cell encode/decode round-trip. Encoding a coordinate to a cell and decoding the cell back to its centre should land inside the same cell. A round-trip that jumps cells signals a resolution or API-version mismatch (a frequent cause when mixing H3 v3 and v4 code paths).

import h3


def assert_h3_roundtrip(lat: float, lon: float, resolution: int) -> None:
    """Assert encode -> decode -> re-encode returns the same H3 cell."""
    cell = h3.latlng_to_cell(lat, lon, resolution)
    c_lat, c_lon = h3.cell_to_latlng(cell)
    recell = h3.latlng_to_cell(c_lat, c_lon, resolution)
    if cell != recell:
        raise AssertionError(
            f"H3 round-trip failed at res {resolution}: {cell} != {recell}"
        )

Cell-coverage sanity. The set of occupied cells should be spatially contiguous with the input footprint and free of stray cells on the far side of the globe—a classic symptom of swapped latitude/longitude arguments. Assert that every occupied cell centre lies within a buffered bounding box of the input.

import h3
import pandas as pd
from shapely.geometry import Point, box


def assert_cells_within_footprint(
    binned: pd.DataFrame, min_lon: float, min_lat: float,
    max_lon: float, max_lat: float, pad_deg: float = 0.05,
) -> None:
    """Assert every occupied H3 cell centre falls inside the padded bbox."""
    if binned.empty:
        raise ValueError("No occupied cells to validate.")
    footprint = box(min_lon - pad_deg, min_lat - pad_deg,
                    max_lon + pad_deg, max_lat + pad_deg)
    for cell in binned["cell"]:
        c_lat, c_lon = h3.cell_to_latlng(cell)
        if not footprint.contains(Point(c_lon, c_lat)):
            raise AssertionError(
                f"Cell {cell} at ({c_lat:.4f}, {c_lon:.4f}) lies outside "
                "the input footprint — check for swapped lat/lon."
            )

Bandwidth units check. Before any kernel density estimate, assert that the coordinates are in a projected metric CRS and that the bandwidth is a sane number of metres. This single guard prevents the degree-space bandwidth bug outright.

import geopandas as gpd


def assert_metric_bandwidth(gdf: gpd.GeoDataFrame, bandwidth_m: float) -> None:
    """Assert coordinates are in a metric CRS and bandwidth is in metres.

    A KDE bandwidth must be a physical distance. Degree-space values
    (typically < 1.0) on a geographic CRS are almost always a mistake.
    """
    if gdf.crs is None:
        raise ValueError("GeoDataFrame has no CRS; cannot validate bandwidth units.")
    if not gdf.crs.is_projected:
        raise ValueError(
            f"CRS {gdf.crs.srs!r} is geographic. Reproject to a metric CRS "
            "(UTM zone or equal-area) before setting a KDE bandwidth."
        )
    if not (1.0 <= bandwidth_m <= 100_000.0):
        raise ValueError(
            f"Bandwidth {bandwidth_m} m is outside a plausible metric range; "
            "did you pass a value in degrees?"
        )

Together these four gates catch the overwhelming majority of aggregation defects: lost points, encode/decode drift, transposed coordinates, and degree-space smoothing. Run them on every batch, not just once during development.

Related

• Discrete Global Grid Systems — H3, S2, and geohash cell indexing for binning pings into density grids
• Map-Matching to Road Networks — snapping noisy traces to road edges with hidden Markov models
• Origin-Destination Flow Matrices — building OD matrices from segmented trip endpoints
• Kernel Density Surfaces — continuous intensity fields with metric-CRS bandwidth discipline
• Movement Pattern Extraction & Trajectory Analysis — the upstream stage that supplies segments and stay points

Back to Home