Aggregating Sensor Data with H3 Cells

This page aggregates a digital twin’s sensor readings into H3 hexagonal cells — choosing a resolution from the sensor spacing rather than by habit, binning and aggregating per cell and time bucket, coarsening through the hierarchy for zoomed-out views, and rendering the result as an overlay on the tiled city, with source coordinates in EPSG:25832 converted to WGS84 for indexing.

Why you hit this

A twin with ten thousand sensors — traffic counters, air quality, noise, occupancy, temperature — cannot show them as ten thousand markers, and their raw positions are not what anybody asks about. The questions are aggregate and spatial: which streets are loudest between seven and nine, where does particulate matter exceed a threshold, how does occupancy vary by district. Binning into a fixed cell system answers all of them with one aggregation, and H3’s hexagons have the property that matters for this: every neighbour is the same distance away, so a diffusion, a gradient or a neighbourhood average behaves consistently in all directions. The alternatives are compared in choosing between S2, H3 and geohash for 3D data.

Prerequisites

  • Python 3.10+ with h3>=4.0, pandas>=2.0, geopandas>=0.14, pyproj>=3.6, shapely>=2.0.
  • Sensor readings with a position, a timestamp and a value. Positions in a projected CRS need converting to WGS84, because H3 indexes latitude and longitude.
  • A sense of the sensor spacing and of the smallest area a reading should be attributed to — both feed the resolution choice.

Step-by-Step

1. Choose the resolution from the sensor spacing

python
import h3

def resolution_table(resolutions=range(7, 14)):
    rows = []
    for r in resolutions:
        area_m2 = h3.average_hexagon_area(r, unit="m^2")
        edge_m = h3.average_hexagon_edge_length(r, unit="m")
        rows.append({"resolution": r,
                     "avg_area_m2": round(area_m2, 1),
                     "avg_edge_m": round(edge_m, 1),
                     "cells_per_km2": round(1e6 / area_m2, 1)})
    return rows

for row in resolution_table():
    print(f"res {row['resolution']:>2}  area {row['avg_area_m2']:>12,.1f} m²  "
          f"edge {row['avg_edge_m']:>8.1f} m  {row['cells_per_km2']:>8.1f} cells/km²")

The rule that works is to pick the resolution whose edge length is close to the sensor spacing. Cells much smaller than the spacing are mostly empty, so the map becomes a scatter of isolated hexagons; cells much larger average away the variation the sensors were installed to measure. For traffic counters every 200–400 m along streets, resolution 9 or 10 is right; for air quality stations several kilometres apart, resolution 7; for indoor occupancy sensors, resolution 13 or finer.

2. Index the readings

python
import pandas as pd
from pyproj import Transformer

to_wgs84 = Transformer.from_crs(25832, 4326, always_xy=True)
RES = 10

def index_readings(df, res=RES):
    lon, lat = to_wgs84.transform(df["easting"].to_numpy(), df["northing"].to_numpy())
    df = df.assign(lon=lon, lat=lat)
    df["cell"] = [h3.latlng_to_cell(la, lo, res) for la, lo in zip(lat, lon)]
    return df

readings = pd.read_parquet("data/noise_readings_2026-09.parquet")
readings = index_readings(readings)
print(f"{len(readings):,} readings, {readings['cell'].nunique():,} distinct cells at res {RES}")
print(readings.head(3)[["sensor_id", "timestamp", "value_db", "cell"]])

H3 takes latitude first and longitude second, which is the opposite order from most projected-coordinate code and the source of a good share of H3 bugs: swap them and every reading lands in the Indian Ocean or the Sahara, both of which are conveniently far from any plausible city and therefore obvious. Less obvious is a partial swap in one code path, which puts a subset of readings somewhere implausible while the rest are fine.

3. Aggregate per cell and time bucket

python
def aggregate(df, freq="1h", value="value_db"):
    df = df.copy()
    df["bucket"] = pd.to_datetime(df["timestamp"]).dt.floor(freq)
    grouped = df.groupby(["cell", "bucket"]).agg(
        n=(value, "size"),
        mean=(value, "mean"),
        p95=(value, lambda s: s.quantile(0.95)),
        max=(value, "max"),
        sensors=("sensor_id", "nunique"),
    ).reset_index()
    grouped["mean"] = grouped["mean"].round(2)
    grouped["p95"] = grouped["p95"].round(2)
    return grouped

agg = aggregate(readings)
print(agg.sort_values("p95", ascending=False).head(5))
print(f"{len(agg):,} cell-hour rows")

Keeping the sensor count per cell alongside the statistics is what makes the result honest. A cell with one sensor and a mean of 71 dB reports that sensor; a cell with six sensors and the same mean reports a neighbourhood. Rendering both identically invites a reader to over-interpret the first, and a map that greys out single-sensor cells at coarse resolutions is usually the right presentation.

4. Coarsen carefully: the hierarchy does not nest exactly

python
def coarsen(agg, from_res=RES, to_res=8):
    a = agg.copy()
    a["parent"] = [h3.cell_to_parent(c, to_res) for c in a["cell"]]
    out = a.groupby(["parent", "bucket"]).agg(
        n=("n", "sum"),
        mean=("mean", "mean"),               # unweighted: see the caveat below
        weighted_mean=("mean", lambda s: None),
        p95=("p95", "max"),
        child_cells=("cell", "nunique"),
    ).reset_index()
    # weight the mean by reading count, which the naive mean of means does not
    wm = (a.assign(prod=a["mean"] * a["n"]).groupby(["parent", "bucket"])
            .agg(num=("prod", "sum"), den=("n", "sum")).reset_index())
    out = out.merge(wm, on=["parent", "bucket"])
    out["weighted_mean"] = (out["num"] / out["den"]).round(2)
    return out.drop(columns=["num", "den"])

coarse = coarsen(agg)
print(coarse.sort_values("n", ascending=False).head(4))

Two traps live in that function. The first is the mean of means, which is wrong whenever cells have different reading counts — the weighted version is one extra group-by and is the number to publish. The second is H3’s hierarchy itself: hexagons cannot tile a hexagon, so a resolution-9 cell is not exactly the union of seven resolution-10 cells. cell_to_parent gives the cell whose centre contains the child’s centre, and a child straddling a parent boundary is assigned wholly to one parent. For aggregation of point readings that is harmless; for aggregating areas — a land-cover share, a footprint total — it introduces an error of a few percent per level, and the right approach there is to re-aggregate from the source rather than to roll up.

Why H3 parents and children do not nest exactly A coarse hexagon overlaid with the seven finer hexagons whose centres fall inside it. The finer cells cover most of the parent and spill over its edges, and slivers of the parent are covered by cells assigned to neighbouring parents. Point readings are unaffected because each point has one cell at each resolution, while area quantities acquire a few percent of error per level. a parent hexagon and the seven children whose centres it contains children spill outside the parent; slivers belong to neighbours points: exact at every resolution areas: a few percent per level Roll up point counts freely; re-aggregate areas from the source.
The hierarchy is approximate by construction, which is a property to plan around rather than a bug to work around.
From readings to an overlay Readings with projected coordinates are transformed to longitude and latitude, indexed to cells at the working resolution, aggregated per cell and hour with a sensor count, optionally coarsened for zoomed-out views, converted to hexagon polygons and rendered as an extruded overlay on the tiled city. Area quantities are re-aggregated from the readings rather than rolled up. readingsE, N, value to lat/lon,index at res 10 aggregate percell and hour coarsen forzoomed views hexagonoverlay area quantities: re-aggregate from readings, never roll up through the hierarchy
One indexing pass serves every later view; the only branch is whether a quantity may be rolled up or must be recomputed.

5. Turn cells into geometry for the overlay

python
import geopandas as gpd
from shapely.geometry import Polygon

def cells_to_gdf(rows, cell_col="cell", crs=4326):
    geoms, records = [], []
    for row in rows.itertuples(index=False):
        cell = getattr(row, cell_col)
        boundary = h3.cell_to_boundary(cell)             # [(lat, lon), …]
        geoms.append(Polygon([(lon, lat) for lat, lon in boundary]))
        records.append(row._asdict())
    gdf = gpd.GeoDataFrame(records, geometry=geoms, crs=crs)
    gdf["area_m2"] = gdf.to_crs(25832).geometry.area.round(1)
    return gdf

peak = agg[agg["bucket"] == agg["bucket"].max()]
hexes = cells_to_gdf(peak)
hexes.to_file("web/noise_hex_res10.geojson", driver="GeoJSON")
print(hexes[["cell", "n", "mean", "p95", "area_m2"]].head(4))
print(f"area range {hexes['area_m2'].min():,.0f}–{hexes['area_m2'].max():,.0f} m²")

cell_to_boundary returns latitude-longitude pairs, so building a Shapely polygon needs them reversed — the same order trap as indexing, in the opposite direction. The area range printed at the end is worth looking at once: H3 cells vary in area with latitude and position within the icosahedral face, typically by a few percent within a city and much more across a country. Any per-area quantity — a density, a rate per hectare — has to divide by each cell’s own area rather than by the resolution’s average.

6. Render it over the tiled city

javascript
const hexes = await Cesium.GeoJsonDataSource.load("/web/noise_hex_res10.geojson", {
  clampToGround: false,
});
viewer.dataSources.add(hexes);

const ramp = [[55, "#1f6b8a"], [60, "#4f7a4d"], [65, "#c46a3d"], [70, "#b0413e"]];
for (const entity of hexes.entities.values) {
  const p95 = entity.properties.p95.getValue();
  const stop = ramp.findLast(([threshold]) => p95 >= threshold) ?? ramp[0];
  entity.polygon.material = Cesium.Color.fromCssColorString(stop[1]).withAlpha(0.55);
  entity.polygon.extrudedHeight = 4 + (p95 - 50) * 2.5;      // exaggerate for readability
  entity.polygon.height = 0;
  entity.polygon.outline = false;
}

Extruding the hexagons by the value turns the overlay into a readable surface over the city rather than a flat choropleth competing with the buildings for attention. Stating the exaggeration in the legend matters, because a 60 m tall hexagon representing 68 dB will otherwise be read as a building.

For a layer that changes continuously — live sensor feeds rather than an hourly aggregate — the geometry should be built once and only the values updated, as in streaming live sensor updates onto tilesets.

Matching resolution to sensor spacing Three panels of the same street network. At resolution 8, with cells about 1.2 kilometres across, all sensors fall into two cells and the variation is averaged away. At resolution 10, with cells about 130 metres across, each cell holds one to three sensors and street-level variation is visible. At resolution 13, with cells nine metres across, almost every cell is empty and the map is a scatter of dots. res 8: one cell, no variation res 10: street-level variation res 13: mostly empty The useful resolution has an edge length close to the sensor spacing.
Resolution is not a quality setting: too fine is as unreadable as too coarse, and the sensor spacing decides which is which.

Expected Output & Verification

text
res  7  area    5,161,293.4 m²  edge  1406.5 m       0.2 cells/km²
res  8  area      737,327.6 m²  edge   531.4 m       1.4 cells/km²
res  9  area      105,332.5 m²  edge   200.8 m       9.5 cells/km²
res 10  area       15,047.5 m²  edge    75.9 m      66.5 cells/km²
res 11  area        2,149.6 m²  edge    28.7 m     465.2 cells/km²
res 12  area          307.1 m²  edge    10.8 m   3,256.2 cells/km²
res 13  area           43.9 m²  edge     4.1 m  22,793.6 cells/km²
1,204,882 readings, 3,118 distinct cells at res 10
        cell              bucket   n   mean    p95   max  sensors
0  8a1fb46622dffff 2026-09-15 08:00  84  71.42  78.10  81.2        3
3,118 distinct cells at res 10
104,204 cell-hour rows
area range 14,802–15,244 m²

Verify the indexing before trusting any aggregate, because a coordinate-order error produces a perfectly consistent map of the wrong place:

python
def verify_indexing(df, res=RES, sample=2000, tolerance_m=200.0):
    from pyproj import Geod
    geod = Geod(ellps="WGS84")
    s = df.sample(min(sample, len(df)), random_state=7)
    dists = []
    for row in s.itertuples(index=False):
        lat_c, lon_c = h3.cell_to_latlng(row.cell)
        _, _, d = geod.inv(row.lon, row.lat, lon_c, lat_c)
        dists.append(d)
    dists = pd.Series(dists)
    return {"checked": len(dists), "max_offset_m": round(float(dists.max()), 1),
             "median_offset_m": round(float(dists.median()), 1),
             "within_cell": bool(dists.max() < tolerance_m)}

check = verify_indexing(readings)
print(check)
assert check["within_cell"], "readings are not inside their cells: latitude/longitude swapped?"

Every reading must lie within its own cell, so the distance from the reading to its cell centre cannot exceed the cell’s circumradius — about 76 m at resolution 10. A median offset of thousands of kilometres is the swapped-order bug; an offset of a few hundred metres on some readings means those were indexed at a different resolution.

Then verify the aggregation preserves the data:

python
assert int(agg["n"].sum()) == len(readings), "readings lost in aggregation"
recomputed = (readings.assign(bucket=pd.to_datetime(readings["timestamp"]).dt.floor("1h"))
              .groupby(["cell", "bucket"])["value_db"].mean().round(2))
merged = agg.set_index(["cell", "bucket"])["mean"]
assert merged.sub(recomputed).abs().max() < 0.011, "aggregated means do not reproduce"
print("aggregation reproduces the source readings")

Performance Notes

  • Indexing is about a microsecond per point in the C-backed library, so a million readings index in a second or two. Vectorise by list comprehension rather than apply, which is several times slower.
  • Aggregate in the database when the data lives there. DuckDB and PostGIS both have H3 extensions, and grouping ten million rows in SQL beats moving them into pandas.
  • Store the cell as a string or a 64-bit integer, not both. The integer form is half the size and sorts usefully; the string form is what most APIs expect. Pick one for storage and convert at the edges.
  • Pre-compute the boundary geometry per cell once and cache it: it is fixed for a given cell, and re-deriving it per render is wasted work.
  • Coarsen from the source for area quantities, from the aggregate for counts and means. The first is correct and slower; the second is fast and approximate.

Common Errors

Every reading lands in the ocean. Latitude and longitude were passed in the wrong order. H3 takes latitude first.

Cell areas differ from the resolution’s average. Expected: H3 cells vary in area. Divide by each cell’s own area for any per-area statistic.

A coarse view shows higher values than any fine cell. The coarse layer used max of the children’s p95, which is a legitimate choice and not the same statistic as the parent’s own p95. Name the statistic in the legend.

Twelve pentagons behave oddly. H3 has twelve pentagonal cells per resolution, inherited from the icosahedron, and they have five neighbours instead of six. Any neighbourhood computation has to tolerate that; grid_disk does it correctly, and hand-rolled neighbour arithmetic does not.

Counts do not add up after a roll-up. Cells were rolled up twice, or a cell straddling a parent boundary was double-counted by a manual containment test. Use cell_to_parent, which assigns each child to exactly one parent.

Frequently Asked Questions

H3 or a square grid?

H3 for anything involving neighbourhoods, flow or diffusion, because all six neighbours are equidistant. A square grid — or the quadkey scheme the tiling already uses — for anything that has to align with tiles, imagery or existing raster products. Mixing them in one twin is normal: hexagons for analysis, quadkeys for delivery.

Can H3 index three dimensions?

No. It is a surface tessellation, so a sensor at street level and one on a roof in the same hexagon share a cell. Where the vertical matters, carry the height as an attribute and bin it separately, or use a 3D scheme as discussed in octree indexing point clouds with Morton codes.

How do I join hexagon aggregates to buildings?

By spatial join, not by cell arithmetic: a building can overlap several cells, and attributing it to the cell containing its centroid is a choice worth stating. For per-building values, aggregate to buildings directly and use hexagons only for the area view.

Back to Spatial Indexing and Tiling Schemes for 3D Data.