Discrete Global Grid Systems

A discrete global grid system (DGGS) tessellates the Earth into a hierarchy of addressable cells so that continuous movement coordinates collapse into a finite set of integer keys you can group, join, and aggregate in O(1).

Prerequisites Checklist

This topic sits inside Spatial Indexing & Density Aggregation and assumes your points have already passed the cleaning stage. Confirm the following before you encode a single cell:

• Python ≥ 3.10 with pandas ≥ 2.0 and geopandas ≥ 0.14 for vectorized DataFrame and geometry operations.
• h3 ≥ 4.0 (the h3-py v4 API — functions are named h3.latlng_to_cell, h3.cell_to_boundary, h3.grid_disk; the v3 names such as geo_to_h3 will raise AttributeError).
• s2sphere or s2geometry if you plan to evaluate S2 cells, and python-geohash for geohash strings.
• Cleaned input. Each point must carry lat (float64), lon (float64), a UTC ts, and a stable device_id. Encoding raw, undeduplicated GPS produces phantom hotspots from multipath scatter — run the point stream through GPS Precision & Error Handling first.
• Two coordinate contexts. Grid encoding consumes WGS84 (EPSG:4326) degrees because H3 and S2 are defined on the sphere. Any distance, dwell, or per-area normalisation must be computed in a metric projected CRS, following Coordinate Reference System Mapping. Keep the two apart in your schema so you never feed metres into a cell encoder.

Grid Distortion Taxonomy

Every grid trades one distortion for another. Choosing a system means choosing which artefact you can tolerate for your analysis. Naming the source up front tells you which correction — area normalisation, resolution change, or a switch of grid family — actually addresses your problem.

Deterministic Pipeline Overview

The five-stage order below is the standard production sequence. Each stage is fully vectorized — no Python-level row iteration — and each preserves the raw point IDs so a cell count can always be traced back to its contributing pings.

1. Choose the grid and resolution. Decide between H3, S2, geohash, and quadkey from the comparison of grid systems for mobility density work, then fix a resolution from the calibration table below. This is the only stage that requires human judgement.
2. Encode cell IDs. Assign every cleaned WGS84 point to its cell with the library’s coordinate-to-cell call. H3 and S2 read spherical lat/lng — never projected metres.
3. Aggregate by cell. Group by cell ID and reduce to ping_count, unique_devices, and a density normalised by the geodesic cell area. This is a single vectorized groupby.
4. Neighbour and k-ring operations. Use h3.grid_disk and h3.grid_ring to build smoothing kernels, contiguity edges for spatial autocorrelation, and adjacency for origin-destination flow matrices.
5. Emit the cell table. Attach each cell’s boundary polygon with h3.cell_to_boundary, and write a cell-indexed GeoParquet table that downstream density and flow stages consume directly.

The complete, runnable version of stages 2–5 for H3 is documented in aggregating GPS pings into H3 hexagon grids; the walkthrough below shows the core encode-and-aggregate step.

Implementation Walkthrough

The function below assigns H3 cells to a cleaned GeoDataFrame and aggregates ping counts and unique devices per cell. H3 encoding is done in WGS84 lat/lng because the grid is defined on the sphere; the density normalisation uses h3.cell_area, which returns a geodesic area in km², so it is metric without ever touching a projected CRS. The one place a projected metric CRS would enter — a within-cell dwell or travel-distance metric — is flagged in the code.

from __future__ import annotations

import logging

import geopandas as gpd
import h3
import numpy as np
import pandas as pd

logger = logging.getLogger(__name__)


def aggregate_points_to_h3(
    gdf: gpd.GeoDataFrame,
    resolution: int = 9,
    device_col: str = "device_id",
    lat_col: str = "lat",
    lon_col: str = "lon",
) -> pd.DataFrame:
    """
    Assign H3 cells to a cleaned point GeoDataFrame and aggregate per cell.

    Parameters
    ----------
    gdf : GeoDataFrame with columns lat (float64), lon (float64), and device_col.
          Coordinates MUST be WGS84 (EPSG:4326) — H3 is defined on the sphere and
          encoding from projected metres yields silently wrong cells.
    resolution : H3 resolution 0–15. See the calibration table for mode guidance.
    device_col : Column identifying the unique moving entity.
    lat_col, lon_col : Column names for latitude and longitude in degrees.

    Returns
    -------
    DataFrame indexed by h3 cell with columns:
      ping_count      — number of points in the cell
      unique_devices  — distinct device_col values in the cell
      area_km2        — geodesic cell area (km²) from h3.cell_area
      density_per_km2 — ping_count / area_km2

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

    if gdf.empty:
        logger.warning("aggregate_points_to_h3 received an empty frame; returning empty.")
        return pd.DataFrame(
            columns=["ping_count", "unique_devices", "area_km2", "density_per_km2"]
        )

    # ── Ensure WGS84 for H3 encoding ─────────────────────────────────────────
    # H3 reads spherical lat/lng. If the frame arrived in a metric CRS from the
    # cleaning stage, reproject a lat/lon view back to EPSG:4326 for encoding.
    # Distances/dwell (not computed here) would instead use the metric CRS.
    if gdf.crs is not None and gdf.crs.to_epsg() != 4326:
        wgs = gdf.to_crs(4326)
        lats = wgs.geometry.y.to_numpy()
        lons = wgs.geometry.x.to_numpy()
    else:
        lats = gdf[lat_col].to_numpy(dtype="float64")
        lons = gdf[lon_col].to_numpy(dtype="float64")

    # ── Stage 2: vectorized cell encoding ────────────────────────────────────
    # h3.latlng_to_cell is scalar; a comprehension over numpy arrays stays fast
    # and avoids DataFrame.iterrows entirely.
    cells = np.fromiter(
        (h3.latlng_to_cell(lat, lon, resolution) for lat, lon in zip(lats, lons)),
        dtype=object,
        count=len(lats),
    )

    work = pd.DataFrame({"h3": cells, device_col: gdf[device_col].to_numpy()})

    # ── Stage 3: aggregate by cell ───────────────────────────────────────────
    agg = work.groupby("h3").agg(
        ping_count=(device_col, "size"),
        unique_devices=(device_col, "nunique"),
    )

    # Geodesic cell area (km²) — metric without a projected CRS.
    agg["area_km2"] = agg.index.map(lambda c: h3.cell_area(c, unit="km^2"))
    agg["density_per_km2"] = agg["ping_count"] / agg["area_km2"]

    logger.info(
        "Aggregated %d points into %d H3 cells at resolution %d",
        len(work), len(agg), resolution,
    )
    return agg.sort_values("ping_count", ascending=False)

The output is a cell-indexed table you can join to any other cell-keyed dataset in constant time, or hand straight to the neighbour operations that build flow adjacency and smoothing kernels.

Geometric and Mathematical Grounding

Hexagon versus square adjacency. On a square grid, a cell has four neighbours across its edges and four across its corners. The edge neighbours sit one cell-width away; the corner neighbours sit √2 ≈ 1.41 cell-widths away. Any operation that treats “the eight neighbours” as equivalent — a smoothing kernel, a spreading model, a contiguity weight — therefore over-weights the diagonal directions and produces a subtly star-shaped result. A hexagon has exactly six neighbours, each sharing a full edge and each at the same centre-to-centre distance. The first ring is a true isotropic distance band, which is why hexagonal grids are the default for kernel density and diffusion over movement data.

Aperture-7 refinement. H3 uses an aperture-7 hierarchy: each parent hexagon is subdivided so that a child resolution has roughly seven times as many cells, and therefore about one-seventh of the area, of its parent. The refinement is not a clean tiling — child hexagons are rotated relative to the parent, so a child can straddle a parent boundary. This is why roll-ups via cell_to_parent are lineage-accurate but leak roughly 1% of points across parent edges. The average cell area collapses by a factor of seven per level:

Pentagons. Twelve pentagon cells exist at every resolution — an unavoidable consequence of closing an icosahedron-based grid over a sphere. Pentagons have five neighbours rather than six, so any k-ring code that hard-codes “6” will miscount near them. In practice the pentagons sit over ocean and remote points by design, but production code should still call h3.is_pentagon before assuming hexagonal adjacency in a study area near one.

Calibration & Parameter Tuning

Resolution is the single most consequential parameter. Pick it by matching the cell edge to both your positional error floor and the scale of the phenomenon — a cell much smaller than your GPS noise scatters one true location across neighbours, while a cell much larger blurs distinct hotspots together.

Two rules keep the resolution honest. First, if the median inter-cell hop for a single stationary device exceeds zero — that is, a parked vehicle keeps changing cells — the resolution is finer than your noise and you should coarsen it or count unique devices over a time window. Second, if more than ~70% of occupied cells hold a single ping, the grid is too fine for the data volume and any density surface will be dominated by sampling noise. Both diagnostics are cheap to compute during the aggregation pass.

Integration & Compatibility

The cell-indexed table this topic produces is a hub that several downstream stages consume:

• Origin-Destination Flow Matrices: The cell IDs become the row and column keys of the OD matrix. Encoding trip endpoints to the same resolution used here lets you build the matrix as a single groupby over (origin_cell, destination_cell) pairs, and grid_disk supplies the adjacency for spatial smoothing of sparse flows.
• Kernel Density Surfaces: Per-cell counts are the discrete input to a hexagonal kernel smoother. Because H3’s first ring is isotropic, a grid_disk-weighted smooth avoids the diagonal bias that a raster kernel over square pixels would introduce.
• GPS Precision & Error Handling: This is the upstream producer. The cleaned, deduplicated, metric-CRS-validated points it emits are the only safe input to cell encoding — aggregating raw fixes bakes multipath scatter into the grid.
• Coordinate Reference System Mapping: Keep the WGS84 encoding context and the metric analysis context explicitly separated in your schema, as covered there, so an encoder never receives projected metres.

FAQ

Should I encode H3 cells from projected metres or from WGS84 lat/lng?

Encode from WGS84 lat/lng. H3, S2, and geohash are all defined on geographic coordinates, so their coordinate-to-cell functions expect degrees. Feeding projected UTM metres into h3.latlng_to_cell produces silently wrong cells. Reserve the metric CRS for distance, dwell, and area work — cell areas for density come from h3.cell_area, which is already geodesic.

Why do my geohash aggregations look distorted at high latitude?

Geohash and quadkey cells are rectangles in a lat/lng or Mercator space, so their ground area shrinks toward the poles. A geohash-6 cell that covers roughly 1.2 km × 0.6 km near the equator collapses to under half that east-west span at 60° latitude, so counts per cell are not comparable across latitudes without explicit area normalisation. H3 keeps a near-uniform cell area worldwide, which is why it is preferred for density work.

What H3 resolution should I use for urban mobility density?

Resolution 9 (≈0.10 km², ~174 m edge) is the workhorse for city-scale pedestrian and micro-mobility density. Drop to resolution 8 (≈0.74 km²) for vehicle origin-destination zones, and go to resolution 10–11 for kerbside analysis. Match the cell edge to your positional error: cells much smaller than the GPS noise floor scatter one true location across several neighbouring hexagons.

How do hexagons help compared with square grids?

Every hexagon has six neighbours all sharing an edge and all equidistant from the centre, so a k-ring is an isotropic distance band. A square cell has four edge neighbours plus four diagonal neighbours at 1.41× the distance, which biases radial smoothing and flow adjacency along the diagonals. For kernel density and spreading processes over movement data, the hexagon’s uniform adjacency removes a directional artefact that squares introduce.

Can I mix resolutions in one aggregation?

Only through the hierarchy, never by concatenating cell IDs from different resolutions in one column. H3 and S2 both let you roll a fine cell up to its parent with cell_to_parent, so aggregate at the finest resolution you need and derive coarser summaries from it. Comparing raw counts across mixed resolutions is meaningless because the cells cover different areas.

Related

• H3 vs S2 vs Quadkey for Mobility Density Grids — decision guide comparing cell shape, area uniformity, hierarchy, and neighbour cost.
• Aggregating GPS Pings into H3 Hexagon Grids — the complete count / unique-device / dwell pipeline with mass-conservation validation.
• Origin-Destination Flow Matrices — the cell IDs from here become OD matrix keys.
• Kernel Density Surfaces — per-cell counts feed a hexagonal kernel smoother.
• GPS Precision & Error Handling — the upstream cleaning stage that produces safe input points.

Back to Spatial Indexing & Density Aggregation