Aggregating Density to H3 Hexagons

A density surface built from parcel centroids and a smoothing kernel looks good and cannot be defended: the value at any point depends on a bandwidth nobody chose for a regulatory reason. Aggregating to a fixed hexagonal grid gives up the smooth appearance and gains a property that matters more — every cell’s value is a count of things that are actually in it, reproducible from the same inputs by anyone. This guide aggregates dwelling units to H3 cells with area-weighted assignment, chooses the resolution from the parcel fabric rather than by eye, and shows why the hierarchical index makes the joins cheap. It is the indexed alternative in automated density calculation grids.

Prerequisites

Step-by-step

Step 1: Choose the resolution from the parcel size

H3 resolutions step by roughly a factor of seven in area. Choosing one is a trade-off between cells that are smaller than a parcel — which forces the apportionment question — and cells so large that the surface says nothing. The rule that holds up is to pick the finest resolution whose cell area comfortably exceeds the median parcel area, so that most parcels sit inside a single cell.

import h3

def choose_resolution(median_parcel_sqm: float, factor: float = 8.0) -> int:
    """Finest resolution whose average cell is at least `factor` x the median parcel."""
    target = median_parcel_sqm * factor
    for res in range(15, -1, -1):
        if h3.average_hexagon_area(res, unit="m^2") >= target:
            return res
    return 0

For a suburban fabric with a median parcel around 800 square metres, this lands on resolution 10 — average cell area about 15,000 square metres, roughly 130 metres across. A rural fabric of large parcels lands on 9 or coarser; a dense urban fabric of small lots on 11. Record the resolution with the output; a density surface without its resolution is uninterpretable.

Choosing the resolution from the parcel, not from the mapAverage H3 cell area by resolution against a typical suburban parcel, showing the finest resolution whose cell comfortably exceeds it.Median suburban parcel800 sq mH3 resolution 112149 sq mH3 resolution 10 — chosen15047 sq mH3 resolution 9105332 sq mThe selector takes the finest resolution whose average cell is at least eight times the median parcel, so most parcels sit inside a single cell.
Below the parcel size the surface shows parcel boundaries rather than density; far above it, the surface says nothing at all.

Step 2: Index the parcels

H3 works in WGS84, so the indexing step uses geographic coordinates while every area calculation stays in the projected frame. Keeping both frames explicit avoids the usual confusion.

import geopandas as gpd

def index_parcels(parcels_proj, res: int):
    """Cell containing each parcel's representative point, from WGS84 coordinates."""
    pts = parcels_proj.geometry.representative_point().to_crs(4326)
    cells = [h3.latlng_to_cell(p.y, p.x, res) for p in pts]
    return parcels_proj.assign(h3_cell=cells)

representative_point rather than centroid is deliberate: a centroid can fall outside a concave parcel, which places its dwelling units in a cell the parcel does not occupy. The representative point is guaranteed to be inside.

Step 3: Decide whether to split parcels across cells, and be explicit

At a resolution eight times the median parcel, most parcels sit in one cell — but the large ones do not, and large parcels are exactly where the dwelling units are. Two treatments are defensible and they answer different questions.

def aggregate_whole(indexed, unit_col="dwelling_units"):
    """Every parcel's units go to the cell containing its representative point."""
    return indexed.groupby("h3_cell")[unit_col].sum()

def aggregate_weighted(parcels_proj, res, unit_col="dwelling_units"):
    """Units distributed across cells in proportion to the parcel's area in each."""
    cells = set()
    for geom in parcels_proj.geometry.to_crs(4326):
        cells.update(h3.geo_to_cells(geom, res))
    grid = gpd.GeoDataFrame(
        {"h3_cell": sorted(cells)},
        geometry=[h3.cells_to_h3shape([c]) for c in sorted(cells)],
        crs=4326).to_crs(parcels_proj.crs)

    parts = gpd.overlay(parcels_proj[["parcel_id", unit_col, "geometry"]],
                        grid, how="intersection", keep_geom_type=True)
    parts["share"] = parts.area / parts["parcel_id"].map(
        parcels_proj.set_index("parcel_id").area)
    parts["units_here"] = parts[unit_col] * parts["share"]
    return parts.groupby("h3_cell")["units_here"].sum()

Whole-parcel assignment answers “how many units are administered from within this cell” and preserves integer counts. Area-weighted assignment answers “how many units are physically in this cell” and produces fractions. For a density standard that is applied per parcel, whole-parcel is usually the honest choice; for a surface used to describe an area, weighted is. What is not acceptable is choosing per-run.

Two assignment policies, answering two different questionsWhole-parcel assignment against area-weighted assignment, compared by the question each answers and what it does to the counts.Whole parcelArea weightedQuestion answeredAdministered from this cellPhysically in this cellCountsStay integersBecome fractionsLarge parcelsAll units to one cellSplit across every cell touchedBest forA standard applied per parcelDescribing an area
What is not acceptable is choosing per run — the two produce different surfaces from the same parcels.

Step 4: Convert counts to density, in the cell’s real area

The average area for a resolution is an average — real H3 cells vary by a few percent with latitude and pentagon adjacency, and using the nominal figure introduces an error that varies systematically across a county.

SQM_PER_ACRE = 4046.8564224

def cell_density(counts, res: int, crs):
    cells = counts.index.tolist()
    grid = gpd.GeoDataFrame(
        {"h3_cell": cells, "units": counts.values},
        geometry=[h3.cells_to_h3shape([c]) for c in cells], crs=4326).to_crs(crs)
    grid["area_acres"] = grid.area / SQM_PER_ACRE if crs.axis_info[0].unit_name.startswith(
        "metre") else grid.area / 43560.0
    grid["units_per_acre"] = grid["units"] / grid["area_acres"]
    return grid

Computing the area from each cell’s projected geometry costs nothing and removes the question entirely. The unit handling here is the same trap covered in unit conversion pitfalls in setback thresholds, which is why the branch on the axis unit is explicit rather than assumed.

Step 5: Use the hierarchy instead of re-aggregating

The property that makes H3 worth the trouble is that a cell’s parent at a coarser resolution is derivable from its index — no spatial operation required. A coarser surface is a string operation and a group-by.

def roll_up(counts, from_res: int, to_res: int):
    """Coarser aggregation without touching geometry — the index carries the hierarchy."""
    assert to_res < from_res
    parent = {c: h3.cell_to_parent(c, to_res) for c in counts.index}
    return counts.groupby(counts.index.map(parent)).sum()

Producing resolutions 10, 9 and 8 from one indexed pass takes milliseconds, where the equivalent with square grids or administrative polygons is three spatial joins. That is also why cross-jurisdiction comparison works: two counties indexed at the same resolution join on the cell identifier directly, with no geometry involved at all.

Why the index is worth the troubleA cell parent at a coarser resolution is derivable from the index itself, so a coarser surface is a string operation and a group-by rather than a spatial join.Index each parcel oncerepresentative point, never the centroidAggregate at the working resolutiona group-by, not a spatial joinRoll up by taking parentsthree resolutions in millisecondsJoin across jurisdictions on the idno geometry, no overlay, no tolerance
Two counties indexed at the same resolution join on the cell identifier directly, with no geometry involved at all.

Verification

The aggregation is correct when nothing is lost, nothing is duplicated, and the surface reconciles against a known total.

total_units = parcels["dwelling_units"].sum()
assert abs(counts.sum() - total_units) < 1e-6, "units lost or duplicated in aggregation"

assert counts.index.is_unique, "duplicate cell ids in the aggregation"
assert all(h3.get_resolution(c) == RES for c in counts.index), "mixed resolutions"

# Rolling up must preserve the total too.
assert abs(roll_up(counts, RES, RES - 2).sum() - total_units) < 1e-6

# Spot-check the densest cell against the parcels inside it.
top = grid.sort_values("units_per_acre", ascending=False).iloc[0]
inside = parcels[parcels["h3_cell"] == top["h3_cell"]]
print(f"{top['h3_cell']}: {top['units_per_acre']:.1f} u/ac from "
      f"{len(inside)} parcels totalling {inside['dwelling_units'].sum()} units")

The conservation assertion is the one that matters. Area-weighted aggregation loses units silently whenever a parcel extends past the indexed cell set, and the total is what reveals it.

Common Pitfalls

  • Using centroids on concave parcels. A centroid outside its own parcel assigns units to the wrong cell. Use a representative point.
  • Dividing by the nominal cell area. Real cells vary by a few percent; computing area from the projected geometry costs nothing and removes the error.
  • Mixing assignment policies between runs. Whole-parcel and area-weighted give different surfaces. Pick one per output and record it.
  • Treating the surface as a regulatory standard. A cell density is a description, not a determination — density standards apply per parcel, and a cell that straddles two districts mixes them.
  • Indexing in a projected CRS. H3 expects WGS84 latitude and longitude; feeding it state plane coordinates produces cells in the wrong hemisphere, usually without an error.
  • Forgetting the pentagons. Twelve cells per resolution are pentagons with slightly different areas and neighbour counts. They rarely fall on inhabited land, but code that assumes six neighbours should not assume it silently.

Frequently Asked Questions

Why hexagons rather than a square grid?

Uniform adjacency, mainly: every hexagon has six neighbours at equal distance, where a square grid has four at one distance and four at another, so any neighbourhood operation on squares is anisotropic. The practical advantages are the hierarchical index and the fact that two datasets indexed at the same resolution join without geometry.

How does this compare with a kernel heatmap?

They answer different questions. A kernel surface is a smoothed estimate whose values depend on a bandwidth; a hexagon surface is an exact count within a stated area. For a map that communicates pattern, the kernel is easier to read — generating density heatmaps from parcel centroids using Rasterio covers that. For anything a person may have to defend, the counted surface is the one to use, because every cell’s value can be traced to the parcels that produced it.

Does the cell index belong in the parcel table?

Yes, and at more than one resolution. Storing the resolution 10 index costs sixteen bytes per parcel and turns every subsequent spatial aggregation into a group-by, including partitioning work for chunking county-scale runs by spatial tile. It also stores well in GeoParquet as a sort key.

What happens at a jurisdiction boundary?

Cells straddle it, which is a feature for regional analysis and a problem for anything jurisdiction-specific. Where the output is per-jurisdiction, either clip the cells or report the boundary cells’ split explicitly — silently attributing a straddling cell to one side is the failure to avoid.

Is resolution 10 always right?

No — it happens to suit a typical suburban fabric. Run the resolution selector against your own median parcel area and check the outcome: if more than a few percent of cells contain a single parcel, the resolution is too fine and the surface is showing parcel boundaries rather than density. If the densest cell contains hundreds of parcels, it is too coarse to show anything.

Can this run incrementally?

Well, in fact. Because a parcel’s cell assignment depends only on that parcel, a refresh only needs to reindex changed parcels and re-sum the affected cells. That makes the aggregation cheap enough to run on every ingest rather than as a periodic batch, which in turn makes the surface a live description of the fabric rather than a snapshot somebody has to remember to refresh.

Part of: Automated density calculation grids