Choosing an H3 Resolution for OSM Point Aggregation Jump to heading#

You have a list of OpenStreetMap points — say every amenity=cafe node in a metro area — and you need to bin them into H3 hexagons, but a resolution that is one step too coarse smears two neighbourhoods into one cell and one step too fine leaves most cells holding a single point or none at all.

Prerequisites Jump to heading#

Confirm each item before running the code; a wrong resolution or an out-of-date h3 build is the usual reason a “density” map ends up either a solid blob or a field of ones.

• h3 ≥ 4.0 installed (pip install "h3>=4.0") — the function names below (latlng_to_cell, cell_to_boundary, average_hexagon_area, grid_disk) are the v4 API and differ from the v3 geo_to_h3 naming.
• Python 3.10+ for the list[tuple[float, float]] and dict typing used here.
• A list of (lon, lat) POI coordinates already extracted from a PBF — see the parent Spatial Index Selection guide for why H3 is the right family for this aggregation in the first place.
• Optional: geopandas ≥ 0.14 and shapely ≥ 2.0 if you want to write hexagon polygons out for mapping.
• The WGS 84 assumption from Coordinate Reference Systems in OSM — H3 consumes latitude and longitude in degrees, so no reprojection is needed.

Conceptual minimum Jump to heading#

H3 is a hierarchy of sixteen resolutions, 0 (coarsest) through 15 (finest), and each step subdivides every parent cell into roughly seven children. That factor of seven is the whole decision: moving one resolution finer divides the average cell area by about seven and multiplies the number of occupied cells accordingly. Choosing a resolution is therefore a balance between three quantities that move together — cell area (how much ground each hexagon covers), cell count (how many hexagons your points spread across), and statistical stability (how many points land in a typical cell). Aggregate at too coarse a resolution and every cell is stable but spatially meaningless; aggregate too fine and each cell is spatially precise but statistically noisy because its count is zero, one, or two. The parent Spatial Index Selection guide explains why hexagons beat a planar grid for this — their near-uniform area keeps counts comparable across latitudes — but it leaves open which of the sixteen resolutions to actually use.

The average area shrinks geometrically with resolution, and the average edge length shrinks with its square root:

$$\bar{A}(r) = \frac{A_0}{7^{\,r}}, \qquad \bar{e}(r) \approx \frac{e_0}{\sqrt{7}^{\,r}}$$

where $A_0 \approx 4.36 \times 10^{6}\ \text{km}^2$ is the average hexagon area at resolution 0 and $e_0 \approx 1108\ \text{km}$ its average edge length. In practical terms this puts resolution 6 near $36\ \text{km}^2$ per cell (neighbourhood-scale), resolution 8 near $0.74\ \text{km}^2$ (a few city blocks), and resolution 10 near $15{,}000\ \text{m}^2$ (a single block). A useful rule of thumb is to target a resolution where the typical occupied cell holds enough points to be stable — often a dozen or more — while cells still resolve the spatial variation you care about.

Runnable solution Jump to heading#

This module aggregates a list of (lon, lat) POIs into per-cell counts at a chosen resolution, and, to make the choice evidence-based, sweeps a range of resolutions and reports the cell count, occupied-cell count, and the median points-per-occupied-cell so you can see the trade-off before committing. It targets h3>=4.0 and Python 3.10+.

from __future__ import annotations

import logging
from collections import Counter
from statistics import median

import h3

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger("osm.h3_aggregation")


def aggregate_points(points: list[tuple[float, float]], resolution: int) -> Counter[str]:
    """Bin (lon, lat) POIs into H3 cells and count per cell.

    h3.latlng_to_cell takes (lat, lng) in that order — OSM coordinates are
    stored (lon, lat), so we swap on the way in.
    """
    if not 0 <= resolution <= 15:
        raise ValueError(f"H3 resolution must be 0..15, got {resolution}")
    counts: Counter[str] = Counter()
    for lon, lat in points:
        cell = h3.latlng_to_cell(lat, lon, resolution)
        counts[cell] += 1
    return counts


def sweep_resolutions(
    points: list[tuple[float, float]], candidates: range = range(5, 11)
) -> list[dict[str, float]]:
    """Report the density trade-off across candidate resolutions."""
    report: list[dict[str, float]] = []
    for res in candidates:
        counts = aggregate_points(points, res)
        occupied = len(counts)
        per_cell = median(counts.values()) if counts else 0
        area_km2 = h3.average_hexagon_area(res, unit="km^2")
        report.append(
            {
                "resolution": res,
                "occupied_cells": occupied,
                "median_per_cell": per_cell,
                "avg_cell_area_km2": round(area_km2, 4),
            }
        )
        logger.info(
            "res %2d | cells=%5d | median/cell=%5.1f | avg area=%.4f km^2",
            res, occupied, per_cell, area_km2,
        )
    return report


def cell_ring_counts(counts: Counter[str], cell: str) -> int:
    """Sum a cell and its six neighbours to smooth a sparse fine grid."""
    return sum(counts.get(c, 0) for c in h3.grid_disk(cell, 1))


def cell_polygon(cell: str) -> list[tuple[float, float]]:
    """Return the hexagon boundary as (lat, lng) vertices for mapping."""
    return h3.cell_to_boundary(cell)


if __name__ == "__main__":
    # Example: café POIs as (lon, lat). In practice load these from a PBF.
    pois: list[tuple[float, float]] = [
        (13.404, 52.520), (13.405, 52.519), (13.388, 52.517),
        (13.412, 52.523), (13.401, 52.521), (13.377, 52.516),
    ]
    table = sweep_resolutions(pois, range(6, 11))
    best = max(t for t in table if t["median_per_cell"] >= 10) \
        if any(t["median_per_cell"] >= 10 for t in table) else table[0]
    logger.info("suggested resolution: %d", int(best["resolution"]))

Step-by-step walkthrough Jump to heading#

1. Guard the resolution. aggregate_points rejects anything outside 0–15 up front, because H3 silently accepts only that band and an out-of-range value is a programming error, not data noise.
2. Swap coordinate order. OSM stores (lon, lat), but h3.latlng_to_cell expects (lat, lng); the function swaps on the way in so callers keep the OSM convention everywhere else. Getting this backwards places every café in the wrong hemisphere.
3. Count into a Counter. Each point contributes one to its cell’s tally, giving a {cell_id: count} map — the core aggregation and the input to any density map.
4. Sweep candidate resolutions. sweep_resolutions runs the aggregation at each candidate and records occupied-cell count and the median points per occupied cell, which is the honest signal of statistical stability — the mean is skewed by a few dense cells.
5. Read the average area. h3.average_hexagon_area(res, unit="km^2") reports the ground area per cell at each resolution, letting you connect the abstract number to the real footprint you want to resolve.
6. Smooth a sparse grid when needed. cell_ring_counts sums a cell with its six grid_disk neighbours, a cheap way to stabilize counts at a fine resolution without dropping to a coarser one.
7. Emit polygons for mapping. cell_to_boundary turns a cell id into hexagon vertices you can hand to GeoPandas or a web map to render the aggregation.
8. Pick from evidence. The __main__ block chooses the finest resolution whose median occupied cell still holds at least ten points, encoding the “stable yet detailed” rule as code rather than a guess.

Verification Jump to heading#

Confirm the aggregation is sound before trusting a density map built on it:

• Counts conserve. sum(counts.values()) must equal len(points) — every POI lands in exactly one cell, so a mismatch means points were dropped or double-counted.
• The sweep is monotonic. In the log, occupied-cell count must rise and median-per-cell must fall as resolution increases; if it does not, the input coordinates are likely swapped or clustered on one point.
• Area matches expectation. h3.average_hexagon_area(8, unit="km^2") should print roughly 0.7373; a wildly different number means a stale h3 v3 install or the wrong unit string.
• Cells round-trip. For any cell, h3.cell_to_boundary returns six vertices (twelve global cells are pentagons and return five), and h3.get_resolution(cell) returns the resolution you aggregated at.
• Neighbours resolve. len(h3.grid_disk(cell, 1)) is 7 for a hexagon, confirming the ring smoothing operates on the expected neighbourhood.

Common errors and fixes Jump to heading#

Specification reference Jump to heading#

H3 defines sixteen resolutions (0–15); each finer resolution subdivides a parent cell into approximately seven children, so average cell area decreases by roughly a factor of seven per step. Cell-to-area, cell-to-boundary, and indexing semantics — including the twelve pentagon cells and the latlng_to_cell argument order — are specified in the official H3 documentation at h3geo.org. Consult the Python statistics documentation for the median used to summarize per-cell counts.

Frequently Asked Questions Jump to heading#

Related Jump to heading#

• Spatial Index Selection: R-tree vs H3 vs Quadkey — the parent guide on why and when to choose H3 over an R-tree or a quadkey grid.
• Spatial Indexing for OSM Extracts — building the R-tree you would pair with an H3 aggregation column.
• Coordinate Reference Systems in OSM — why H3 consumes WGS 84 degrees directly with no reprojection.
• Node-Way-Relation Data Model — where the POI node coordinates you aggregate come from.
• Tag Taxonomy & Key-Value Standards — selecting which tagged features (for example amenity=cafe) to aggregate.

Up one level: Spatial Index Selection: R-tree vs H3 vs Quadkey.