H3 vs S2 vs quadkey for mobility density grids

Choosing a grid for mobility density is a one-way door: it fixes your cell shape, your neighbour semantics, and how comparable your counts are across latitudes, and re-binning a billion-row point archive later is expensive. The short version — H3 for density and diffusion, S2 for exact hierarchical covering, quadkey for web-tile interoperability. This page bins the same sample of cleaned points three ways so the trade-offs are concrete rather than theoretical.

Why the choice matters

All three systems belong to the family of grids introduced in Discrete Global Grid Systems, but they differ in the two properties that dominate density work: cell shape (which sets neighbour geometry) and area uniformity (which sets whether raw counts are comparable across space). Get either wrong and your heatmap encodes an artefact of the grid rather than a signal in the movement.

• H3 tiles the globe with hexagons (plus twelve unavoidable pentagons). Every hexagon has six edge-neighbours at equal distance, and cell area is near-uniform worldwide. That isotropy is exactly what a kernel density surface wants.
• S2 projects the sphere onto the six faces of a cube and lays a quadtree on each face. Cells are quadrilaterals with a strict aperture-4 hierarchy and 64-bit integer IDs — excellent for exact covering and database keys, less ideal for isotropic smoothing.
• Quadkey is the Web Mercator tile scheme used by slippy maps. Cells are Mercator squares addressed by a base-4 string. Its ground area shrinks toward the poles as 1/cos(latitude), so it is the least uniform of the three.

Whichever you choose, the points must already be cleaned and metric-CRS-validated per GPS Precision & Error Handling — grid choice cannot rescue noisy input.

Head-to-head comparison

The decisive rows for most density work are neighbour geometry and area uniformity, and on both H3 leads. S2’s strength is hierarchy and integer indexing; quadkey’s is interoperability with tile infrastructure you may already run.

Binning the same points three ways

The script below encodes one identical sample of cleaned WGS84 points into all three grids at comparable resolutions, then reports how many cells each produces and how uniform their areas are. Encoding happens in WGS84 lat/lng for every grid because all three are defined on geographic coordinates; the area comparison uses each library’s own geodesic/analytic area so no projected CRS is needed. Any per-cell distance metric you add later must be computed in a metric CRS — that discipline is unchanged by the grid choice.

from __future__ import annotations

import h3
import numpy as np
import pandas as pd
import s2sphere
import mercantile


def bin_points_three_ways(
    df: pd.DataFrame,
    lat_col: str = "lat",
    lon_col: str = "lon",
    h3_res: int = 9,
    s2_level: int = 13,
    quad_zoom: int = 16,
) -> dict[str, pd.DataFrame]:
    """
    Bin the SAME cleaned WGS84 points into H3, S2, and quadkey cells and
    return one aggregated count table per grid for comparison.

    Parameters
    ----------
    df : DataFrame with lat_col, lon_col in WGS84 degrees (EPSG:4326).
         All three grids encode from spherical lat/lng — never projected metres.
    h3_res : H3 resolution (9 ≈ 0.10 km²).
    s2_level : S2 level (13 ≈ 0.11 km², comparable to H3 res 9).
    quad_zoom : Web Mercator tile zoom (16 ≈ 0.10–0.15 km² near mid-latitude).

    Returns
    -------
    dict mapping grid name → DataFrame indexed by cell id with columns
    ping_count and area_km2.

    Raises
    ------
    ValueError if columns are missing or the frame is empty.
    """
    required = {lat_col, lon_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if df.empty:
        raise ValueError("Input DataFrame is empty.")

    lats = df[lat_col].to_numpy(dtype="float64")
    lons = df[lon_col].to_numpy(dtype="float64")

    # ── H3 encoding (WGS84 → hex cell) ───────────────────────────────────────
    h3_cells = [h3.latlng_to_cell(la, lo, h3_res) for la, lo in zip(lats, lons)]
    h3_tbl = (
        pd.Series(h3_cells, name="cell").value_counts().rename("ping_count").to_frame()
    )
    h3_tbl["area_km2"] = h3_tbl.index.map(lambda c: h3.cell_area(c, unit="km^2"))

    # ── S2 encoding (WGS84 → cube-face quadtree cell) ────────────────────────
    def s2_cell(lat: float, lon: float) -> str:
        ll = s2sphere.LatLng.from_degrees(lat, lon)
        return s2sphere.CellId.from_lat_lng(ll).parent(s2_level).to_token()

    s2_cells = [s2_cell(la, lo) for la, lo in zip(lats, lons)]
    s2_tbl = (
        pd.Series(s2_cells, name="cell").value_counts().rename("ping_count").to_frame()
    )
    # S2 exact cell area in steradians → km² (Earth radius 6371.0088 km).
    R2 = 6371.0088 ** 2
    s2_tbl["area_km2"] = s2_tbl.index.map(
        lambda t: s2sphere.Cell(s2sphere.CellId.from_token(t)).exact_area() * R2
    )

    # ── Quadkey encoding (WGS84 → Web Mercator tile) ─────────────────────────
    def quadkey(lat: float, lon: float) -> str:
        tile = mercantile.tile(lon, lat, quad_zoom)
        return mercantile.quadkey(tile)

    quad_cells = [quadkey(la, lo) for la, lo in zip(lats, lons)]
    quad_tbl = (
        pd.Series(quad_cells, name="cell").value_counts().rename("ping_count").to_frame()
    )
    # Web Mercator tile ground area shrinks with latitude: multiply the
    # equatorial tile area by cos(centre latitude) as an approximation.
    def quad_area_km2(qk: str) -> float:
        tile = mercantile.quadkey_to_tile(qk)
        b = mercantile.bounds(tile)
        lat_c = np.radians((b.north + b.south) / 2.0)
        deg_km = 111.320  # km per degree at the equator
        width = (b.east - b.west) * deg_km * np.cos(lat_c)
        height = (b.north - b.south) * deg_km
        return abs(width * height)

    quad_tbl["area_km2"] = quad_tbl.index.map(quad_area_km2)

    return {"h3": h3_tbl, "s2": s2_tbl, "quadkey": quad_tbl}


def summarise_binning(tables: dict[str, pd.DataFrame]) -> pd.DataFrame:
    """Report cell count, total pings, and area coefficient of variation per grid."""
    rows = []
    for name, tbl in tables.items():
        area = tbl["area_km2"].to_numpy()
        cv = float(np.std(area) / np.mean(area)) if len(area) else float("nan")
        rows.append(
            {
                "grid": name,
                "n_cells": len(tbl),
                "total_pings": int(tbl["ping_count"].sum()),
                "area_cv": round(cv, 4),   # lower = more uniform cell area
                "mean_area_km2": round(float(np.mean(area)), 4),
            }
        )
    return pd.DataFrame(rows).set_index("grid")

Running summarise_binning on a real city sample makes the theory measurable: the area_cv column is near zero for H3, small but non-zero for S2, and largest for quadkey — the quantitative signature of each grid’s area distortion. The total_pings count must be identical across all three grids; if it is not, a point fell outside a grid’s valid domain and needs investigating before you trust any count.

Choose H3 when…

• Your primary output is a density surface, hotspot map, or any diffusion/smoothing process — the isotropic six-neighbour ring is worth more than anything the alternatives offer.
• You want cell counts to be comparable across a wide latitude range without an area correction.
• You want the simplest neighbour API: h3.grid_disk(cell, k) returns a full k-ring in one call. The complete H3 workflow is in aggregating GPS pings into H3 hexagon grids.

Choose S2 when…

• You need an exact, gap-free hierarchical covering of an arbitrary polygon (for example, a city boundary), which S2’s RegionCoverer produces directly.
• You want strict parent-child nesting and compact 64-bit integer keys for a spatial database index, where H3’s approximate aperture-7 nesting would leak points across parents.
• You are already inside an S2-based stack (some ride-hail and geo-database systems standardise on it).

Choose quadkey when…

• You must join mobility density onto existing XYZ / slippy-map raster tiles, or feed a web-map front end that already speaks quadkeys.
• Your analysis stays within a single mid-latitude city, where Mercator area distortion across the study area is small and interoperability outweighs uniformity.
• You value trivial, arithmetic-only cell addressing and tile-pyramid tooling over neighbour isotropy.

Common mistakes and gotchas

• Comparing resolutions across grids by number. H3 res 9, S2 level 13, and quadkey zoom 16 are only roughly area-matched. Always compare by mean cell area (as the script does), not by the raw level integer.
• Ignoring quadkey latitude scaling across cities. A national or continental density map on quadkeys will over-count northern cells unless you divide by the latitude-corrected area.
• Assuming H3 children tile their parent exactly. They do not — the aperture-7 rotation leaks about 1% of points across parent boundaries. If you need exact roll-ups, that is an argument for S2.
• Encoding from projected metres. All three encoders expect WGS84 degrees. Passing UTM coordinates yields cells on the wrong part of the planet with no error raised.

Related

• Discrete Global Grid Systems — the parent overview of the grid family, distortion taxonomy, and encoding pipeline.
• Aggregating GPS Pings into H3 Hexagon Grids — the full H3 count / unique-device / dwell workflow once you have chosen H3.
• Kernel Density Surfaces — where isotropic hexagon neighbours pay off in smoothing.
• GPS Precision & Error Handling — the cleaning stage every grid depends on.

Back to Discrete Global Grid Systems