Aggregating GPS pings into H3 hexagon grids

Turning a raw stream of GPS pings into a hexagon-indexed density table is a four-step job: assign pings to cells → aggregate per cell → build hex polygons → validate mass conservation. Done right, you get a compact GeoDataFrame where every row is one H3 cell carrying a ping count, a unique-device count, and dwell seconds, ready to hand to a heatmap or flow stage. Done carelessly, you get double-counted stationary devices, hexagons drawn on the wrong side of the meridian, and counts that silently drop points on the floor.

Why this needs its own pipeline

Binning looks trivial — “group by cell” — but three failure modes make it worth a dedicated, tested function. First, the encode step is easy to get wrong by feeding projected metres into an encoder that expects degrees. Second, dwell and unique-device aggregation require per-device time ordering that a naïve groupby skips. Third, the polygon construction depends on a vertex order and coordinate convention that h3.cell_to_boundary does not spell out. This deep-dive is the runnable companion to the encoding overview in Discrete Global Grid Systems; if you have not yet decided that H3 is the right grid, start with the comparison of H3, S2, and quadkey.

The input must already be clean. Aggregating pings that still contain GPS drift bakes multipath scatter into your hexagons as phantom neighbouring-cell hits, so run the stream through handling GPS drift in raw trajectory logs before you encode.

Core aggregation pipeline

H3 Ping Aggregation Pipeline Four sequential stages: Assign pings to H3 cells, Aggregate per cell into count, unique devices and dwell, Build hex GeoDataFrame with cell boundaries, and Validate mass conservation, connected by arrows. 1. Assign to H3 cells 2. Aggregate count · devices · dwell 3. Build hex GeoDataFrame 4. Validate mass conservation
  1. Assign pings to H3 cells. Encode each cleaned WGS84 ping to its cell at the chosen resolution. This is a vectorized comprehension over NumPy arrays — never DataFrame.iterrows.
  2. Aggregate per cell. Group by cell for the ping count and the unique-device count; compute dwell by ordering each device’s pings in time first.
  3. Build the hex GeoDataFrame. Attach each occupied cell’s boundary polygon with h3.cell_to_boundary, remembering the (lat, lng) → (lng, lat) swap, and assemble a GeoDataFrame in EPSG:4326.
  4. Validate mass conservation. Assert that the summed cell counts equal the input ping count and that every aggregated device actually appears in the input.

Production-ready Python implementation

The function computes dwell per ping as the time to that device’s next ping, capped so a sparse gap cannot credit a cell with phantom hours. H3 encoding uses WGS84 lat/lng — the grid is spherical, so projected metres would be wrong; if you also derive a travel-distance metric it must be computed in a metric CRS, which this aggregation deliberately does not touch.

PYTHON 
from __future__ import annotations

import logging

import geopandas as gpd
import h3
import numpy as np
import pandas as pd
from shapely.geometry import Polygon

logger = logging.getLogger(__name__)

_EMPTY_COLS = ["h3", "ping_count", "unique_devices", "dwell_s", "area_km2", "geometry"]


def aggregate_pings_to_h3(
    df: pd.DataFrame,
    resolution: int = 9,
    lat_col: str = "lat",
    lon_col: str = "lon",
    ts_col: str = "ts",
    device_col: str = "device_id",
    max_dwell_s: float = 300.0,
) -> gpd.GeoDataFrame:
    """
    Bin cleaned GPS pings into H3 hexagons with count, unique-device, and dwell.

    Parameters
    ----------
    df : DataFrame with lat_col, lon_col (WGS84 degrees, EPSG:4326), a UTC ts_col,
         and device_col. Coordinates must be geographic — H3 encodes from lat/lng,
         not from projected metres.
    resolution : H3 resolution 0–15.
    max_dwell_s : Per-ping dwell is capped here so a long gap between two sparse
                  pings does not credit a cell with phantom dwell time.

    Returns
    -------
    GeoDataFrame (EPSG:4326) indexed 0..n with columns:
      h3, ping_count, unique_devices, dwell_s, area_km2, geometry (hex Polygon).
      An empty input yields an empty GeoDataFrame with the same columns and CRS.

    Raises
    ------
    ValueError if required columns are missing or resolution is out of range.
    """
    required = {lat_col, lon_col, ts_col, device_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if not 0 <= resolution <= 15:
        raise ValueError(f"H3 resolution must be 0–15, got {resolution}")

    # ── Empty-frame guard: return a well-formed empty GeoDataFrame ────────────
    if df.empty:
        logger.warning("aggregate_pings_to_h3 received an empty frame; returning empty.")
        return gpd.GeoDataFrame(
            {c: pd.Series(dtype="float64") for c in _EMPTY_COLS},
            geometry="geometry",
            crs="EPSG:4326",
        )

    work = df[[lat_col, lon_col, ts_col, device_col]].copy()
    work[ts_col] = pd.to_datetime(work[ts_col], utc=True)

    # ── Stage 1: vectorized H3 cell assignment (WGS84 → cell) ─────────────────
    lats = work[lat_col].to_numpy(dtype="float64")
    lons = work[lon_col].to_numpy(dtype="float64")
    work["h3"] = [h3.latlng_to_cell(la, lo, resolution) for la, lo in zip(lats, lons)]

    # ── Dwell: time to each device's next ping, capped, assigned to its cell ──
    work = work.sort_values([device_col, ts_col])
    nxt = work.groupby(device_col)[ts_col].shift(-1)
    dwell = (nxt - work[ts_col]).dt.total_seconds()
    work["dwell_s"] = dwell.clip(lower=0.0, upper=max_dwell_s).fillna(0.0)

    # ── Stage 2: aggregate per cell ──────────────────────────────────────────
    agg = work.groupby("h3").agg(
        ping_count=(device_col, "size"),
        unique_devices=(device_col, "nunique"),
        dwell_s=("dwell_s", "sum"),
    )
    agg["area_km2"] = agg.index.map(lambda c: h3.cell_area(c, unit="km^2"))

    # ── Stage 3: build hex polygons from cell boundaries ─────────────────────
    # h3.cell_to_boundary returns (lat, lng) pairs; Shapely wants (lng, lat).
    def hex_polygon(cell: str) -> Polygon:
        boundary = h3.cell_to_boundary(cell)  # tuple of (lat, lng)
        return Polygon([(lng, lat) for lat, lng in boundary])

    agg = agg.reset_index()
    agg["geometry"] = agg["h3"].map(hex_polygon)

    gdf = gpd.GeoDataFrame(agg, geometry="geometry", crs="EPSG:4326")
    logger.info(
        "Binned %d pings into %d H3 cells at resolution %d",
        len(work), len(gdf), resolution,
    )
    return gdf

Validation block: mass conservation

The single most important check on an aggregation is that it neither invents nor loses observations. Summed cell counts must equal the input ping count, and every device credited to a cell must appear in the raw input. Run this immediately after aggregate_pings_to_h3.

PYTHON 
def validate_h3_aggregation(
    raw: pd.DataFrame,
    gdf: gpd.GeoDataFrame,
    device_col: str = "device_id",
) -> None:
    """
    Assert mass conservation for an H3 aggregation.
    Raises AssertionError on failure; prints a summary on success.
    """
    # 1. Ping mass conservation: no point invented or dropped.
    total_cells = int(gdf["ping_count"].sum()) if len(gdf) else 0
    assert total_cells == len(raw), (
        f"Ping mass not conserved: {len(raw)} in → {total_cells} binned"
    )

    # 2. Unique-device sanity: per-cell distinct devices cannot exceed the
    #    global distinct-device count.
    global_devices = raw[device_col].nunique()
    if len(gdf):
        assert gdf["unique_devices"].max() <= global_devices, (
            "A cell reports more unique devices than exist in the input"
        )

    # 3. Geometry validity: every hexagon must be a closed, valid polygon
    #    with coordinates in WGS84 range (catches the lat/lng swap bug).
    if len(gdf):
        assert gdf.geometry.is_valid.all(), "Invalid hexagon geometry produced"
        bounds = gdf.total_bounds  # minx, miny, maxx, maxy
        assert -180 <= bounds[0] and bounds[2] <= 180, "Longitude out of range — check lat/lng swap"
        assert -90 <= bounds[1] and bounds[3] <= 90, "Latitude out of range — check lat/lng swap"

    # 4. Non-negative dwell.
    if len(gdf):
        assert (gdf["dwell_s"] >= 0).all(), "Negative dwell seconds"

    print(
        f"Validation passed: {len(raw)} pings → {len(gdf)} cells, "
        f"{int(gdf['ping_count'].sum()) if len(gdf) else 0} pings conserved."
    )

Mass conservation is the assertion that catches the widest class of bugs: a silent CRS mismatch that drops points outside the grid, a groupby that discards NaN cells, or a filter applied in the wrong order. If pings in does not equal pings binned, stop and fix the pipeline before trusting any density it produces. These per-cell counts are exactly the input the kernel density surfaces stage smooths, so an error here propagates straight into the heatmap.

Common mistakes and gotchas

  • Feeding projected metres into h3.latlng_to_cell. The encoder expects WGS84 degrees. UTM eastings/northings look like plausible numbers, so no error is raised — the cells simply land in the wrong hemisphere. If your frame arrived in a metric CRS, reproject a lat/lon view back to EPSG:4326 before encoding.
  • Forgetting the (lat, lng) → (lng, lat) swap. h3.cell_to_boundary returns latitude-first pairs, but Shapely’s Polygon wants (x, y) = (lng, lat). Skip the swap and every hexagon is mirrored across the diagonal — the validation block’s bounds check catches it.
  • Double-counting stationary devices. A parked device that jitters across a cell edge between epochs inflates the unique-device count in both cells. Count distinct devices per cell over a time window rather than per raw ping, and consider snapping known dwell locations to a single cell.
  • Uncapped dwell across sparse gaps. If two consecutive pings for a device are ten minutes apart, an uncapped dwell credits the first cell with ten minutes it may never have been occupied. The max_dwell_s cap bounds this; tune it to your sampling interval.
  • Choosing a resolution finer than your noise floor. At resolution 11 on consumer GPS, one true stationary location smears across three or four neighbouring hexagons. Coarsen the resolution or smooth the input first; the resolution guidance in the grid systems overview maps cell edge to transport mode.

Back to Discrete Global Grid Systems

Related in this section