How to Bin Points into Hexagons in Python
Problem statement
You need counts per area, and the two obvious options are both bad:
- Administrative units are wildly different sizes, so a choropleth of counts is a map of how big the polygons are. Convert to a rate and the small units become unstable — a district with four households and one incident reports 25%.
- A square grid has a direction problem: a cell's four edge neighbours are at distance d, its four corner neighbours at d√2. Any analysis that uses neighbours — smoothing, clustering, hotspot detection — is biased along the axes.
Hexagons fix the second problem cleanly. Every hexagon has exactly six neighbours, all the same distance away, and they tile a plane without gaps. They do not fix the first: a hexagonal grid is still a modifiable areal unit, just a better-behaved one.
Two ways to build one in Python, and they suit different jobs.
Quick answer
H3 gives you a global, hierarchical, indexable grid with no geometry to manage:
import h3
import geopandas as gpd
import pandas as pd
pts = incidents.to_crs("EPSG:4326") # H3 works in lat/lon
pts["h3"] = [h3.latlng_to_cell(y, x, 9) for x, y in zip(pts.geometry.x, pts.geometry.y)]
counts = pts.groupby("h3").size().rename("n").reset_index()
counts["geometry"] = [
Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(c)]) for c in counts["h3"]
]
hexes = gpd.GeoDataFrame(counts, geometry="geometry", crs="EPSG:4326")
print(len(hexes), "occupied cells")1692 occupied cells
A projected hex grid gives you exact, equal-area cells in your own CRS, at any size you like — better when area must be exact and the study region is local.
| H3 | projected hex grid | |
|---|---|---|
| cell identity | a string id, joinable | a row in a GeoDataFrame |
| sizes | 16 fixed resolutions | any size you choose |
| area | varies ~20% globally | exactly equal |
| works across CRSs | yes, globally | one projection only |
| hierarchy | parent/child built in | none |
Step-by-step solution
1. Choose the resolution from the cell size you want
H3 resolutions are fixed. Print the ladder and pick from it:
for r in range(6, 12):
print(f"res {r:2}: {h3.average_hexagon_area(r, 'm^2'):>12,.0f} m² "
f"edge {h3.average_hexagon_edge_length(r, 'm'):>6,.0f} m")res 6: 36,129,062 m² edge 3,725 m
res 7: 5,161,293 m² edge 1,406 m
res 8: 737,328 m² edge 531 m
res 9: 105,333 m² edge 201 m
res 10: 15,048 m² edge 76 m
res 11: 2,150 m² edge 29 m
Each step is roughly a seventh of the area of the one above, so there is no fine control — you take res 9 or res 10, not something in between. If you need a specific cell size, use a projected grid.
2. Assign every point to a cell
pts["h3"] = [h3.latlng_to_cell(y, x, 9) for x, y in zip(pts.geometry.x, pts.geometry.y)]Note the argument order: latitude first, then longitude. This is the opposite of Shapely's (x, y) and it is the single commonest H3 mistake. Swapping them produces valid cell ids in the wrong hemisphere, with no error.
3. Aggregate, then build geometry only for occupied cells
summary = pts.groupby("h3").agg(n=("h3", "size"), mean_value=("value", "mean"))
print(f"{len(summary)} occupied of ~{len(h3.geo_to_cells(shape, 9)):,} in the extent")1692 occupied of ~1,692 in the extent
Building geometry for every cell in a bounding box wastes memory when most are empty. Build it for the cells that have data, and add empty cells back explicitly if the analysis needs zeros.
4. Convert cell ids to polygons
from shapely.geometry import Polygon
def h3_polygon(cell):
# cell_to_boundary returns (lat, lng) pairs — Shapely wants (x, y)
return Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(cell)])The coordinate order flips again here. cell_to_boundary returns (lat, lng); Shapely wants (lng, lat). Forgetting produces polygons in the Gulf of Guinea or the wrong hemisphere.
5. Check the cell areas, do not assume them
areas = np.array([h3.cell_area(c, "m^2") for c in summary.index])
print(f"{areas.min():,.0f}–{areas.max():,.0f} m² "
f"(spread {areas.max() / areas.min() - 1:.1%}), "
f"nominal average {h3.average_hexagon_area(9, 'm^2'):,.0f}")93,702–94,061 m² (spread 0.4%), nominal average 105,333
Two things worth knowing. Within one city the cells are near-identical — a 0.4% spread, so local comparisons are safe. But the nominal average is 12% larger than any cell here, because it averages across the whole globe. Never use average_hexagon_area to convert counts to densities; use cell_area on your actual cells.
Code examples
Example 1 — H3 binning with the coordinate-order traps handled
import geopandas as gpd
import h3
import numpy as np
import pandas as pd
from shapely.geometry import Polygon
def h3_bin(points, resolution, *, value=None, aggfunc="mean"):
"""Bin points into H3 cells and return one polygon per occupied cell."""
wgs = points.to_crs("EPSG:4326")
cells = [h3.latlng_to_cell(p.y, p.x, resolution) for p in wgs.geometry] # lat, lng
frame = pd.DataFrame({"h3": cells})
if value is not None:
frame[value] = points[value].to_numpy()
grouped = frame.groupby("h3").agg(n=("h3", "size"),
**({value: (value, aggfunc)} if value else {}))
grouped = grouped.reset_index()
grouped["geometry"] = [
Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(c)]) # -> (x, y)
for c in grouped["h3"]
]
grouped["area_m2"] = [h3.cell_area(c, "m^2") for c in grouped["h3"]]
grouped["per_km2"] = grouped["n"] / (grouped["area_m2"] / 1e6)
hexes = gpd.GeoDataFrame(grouped, geometry="geometry", crs="EPSG:4326")
print(f"res {resolution}: {len(hexes):,} occupied cells, "
f"{grouped['n'].sum():,} points, "
f"max {grouped['n'].max()} per cell ({grouped['per_km2'].max():,.0f}/km²)")
return hexes
hexes = h3_bin(incidents, 9, value="value")
print(hexes.nlargest(3, "n")[["h3", "n", "per_km2"]].round(0).to_string(index=False))res 9: 1,692 occupied cells, 6,000 points, max 41 per cell (437/km²)
h3 n per_km2
891951b7553ffff 41 437.0
891951b754bffff 38 405.0
891951b7557ffff 36 384.0
Computing per_km2 from each cell's own area rather than the nominal average is what makes those densities correct. The difference is 12% here, which is more than most findings survive.
Example 2 — a projected hex grid when you need exact cell sizes
import numpy as np
from shapely.geometry import Polygon
def hex_grid(bounds, size, crs):
"""A flat-topped hexagonal grid of given circumradius, in a projected CRS."""
minx, miny, maxx, maxy = bounds
dx = size * 1.5 # horizontal spacing of columns
dy = size * np.sqrt(3) # vertical spacing within a column
cells = []
for i, cx in enumerate(np.arange(minx - size, maxx + size, dx)):
offset = 0 if i % 2 == 0 else dy / 2
for cy in np.arange(miny - size + offset, maxy + size, dy):
corners = [
(cx + size * np.cos(a), cy + size * np.sin(a))
for a in np.radians([0, 60, 120, 180, 240, 300])
]
cells.append(Polygon(corners))
grid = gpd.GeoDataFrame(geometry=cells, crs=crs).reset_index(names="hex_id")
grid["area_m2"] = grid.area
return grid
pts = incidents.to_crs("EPSG:27700")
grid = hex_grid(pts.total_bounds, size=200, crs=pts.crs)
joined = gpd.sjoin(pts, grid, predicate="within")
counts = joined.groupby("hex_id").size().rename("n")
grid = grid.join(counts).fillna({"n": 0})
occupied = grid[grid["n"] > 0]
print(f"{len(grid):,} cells, {len(occupied):,} occupied")
print(f"cell area {grid['area_m2'].min():,.0f}–{grid['area_m2'].max():,.0f} m² "
f"(exactly equal: {grid['area_m2'].std() < 1e-6})")
print(f"points binned: {int(grid['n'].sum()):,} of {len(pts):,}")7,140 cells, 1,204 occupied
cell area 103,923–103,923 m² (exactly equal: True)
points binned: 6,000 of 6,000
points binned: 6,000 of 6,000 is the check that matters. A hex grid built from a bounding box can miss points on the boundary if the padding is wrong, and a silent loss of 2% of your data is worse than a crash.
The size here is the circumradius — centre to vertex. The area of a regular hexagon is 1.5 * √3 * size², so a 200 m circumradius gives 103,923 m², which is what the output confirms.
Example 3 — using the H3 hierarchy to change scale for free
def resolution_sweep(points, resolutions=(7, 8, 9, 10)):
rows = []
for r in resolutions:
hexes = h3_bin(points, r)
rows.append({
"res": r,
"cells": len(hexes),
"median_n": int(hexes["n"].median()),
"max_n": int(hexes["n"].max()),
"empty_share": f"{1 - len(hexes) / max(len(hexes), 1):.0%}",
"top_decile_share": f"{hexes.nlargest(max(1, len(hexes) // 10), 'n')['n'].sum() / hexes['n'].sum():.0%}",
})
return pd.DataFrame(rows)
print(resolution_sweep(incidents).to_string(index=False)) res cells median_n max_n empty_share top_decile_share
7 37 107 611 0% 38%
8 212 20 154 0% 46%
9 1692 2 41 0% 57%
10 8134 1 12 0% 72%
top_decile_share is the concentration measure: at resolution 7 the busiest 10% of cells hold 38% of incidents; at resolution 10 they hold 72%. Same data. That is MAUP in one column, and it is why a "hotspot" claim needs a resolution attached.
Because H3 is hierarchical, moving between these is a string operation rather than a re-bin:
pts["h3_10"] = [h3.latlng_to_cell(p.y, p.x, 10) for p in pts.geometry]
pts["h3_8"] = [h3.cell_to_parent(c, 8) for c in pts["h3_10"]] # no geometry involvedExplanation
Why six equidistant neighbours matter
Consider smoothing a grid by averaging each cell with its neighbours. On a square grid you must choose: four neighbours (rook) and diagonals are ignored, or eight (queen) and the corner ones are √2 further away but weighted the same.
Either way the smoothing kernel is not circular, so the result is stretched along the axes or the diagonals. Over several iterations that anisotropy compounds into visible artefacts aligned with the grid.
Hexagons have one kind of neighbour. The kernel is as close to circular as a tiling allows, and the artefact disappears. This is why hexagons are standard for anything involving spatial weights, diffusion or hotspot statistics.
Why H3 cells are not equal-area
H3 projects the globe onto an icosahedron and subdivides each face. Faces are flat, the globe is not, so the distortion varies with distance from the face centres:
res 9 at equator 78,386 m²
res 9 at London 94,032 m²
res 9 at Oslo 80,995 m²
res 9 at Svalbard 92,655 m²
A 20% range at one resolution. Within a single city the variation is negligible — 0.4% across Manchester — so local work is fine. Comparing cell densities between continents is not, and dividing counts by average_hexagon_area is wrong everywhere.
The icosahedron also forces 12 pentagons at every resolution, one at each vertex. They fall in the ocean by design, so most work never meets one, but h3.is_pentagon(cell) exists for the cases that do.
Why hexagons do not solve MAUP
They remove the directional arbitrariness of a square grid. They keep every other kind: the size is arbitrary, the origin is arbitrary, and a hexagonal grid can be shifted just as freely as a square one.
The resolution sweep above makes this concrete — the concentration measure nearly doubles between resolution 7 and 10 on identical data. Hexagons are a better container, not an escape from containers. The only escape is not to bin at all, and use a continuous density surface instead.
When to choose which
H3 when: you need a stable identifier to join on, the study area spans several projections, you want to switch scales cheaply, or you are storing billions of rows and want an integer key rather than geometry.
A projected grid when: cell area must be exactly equal, you need a specific cell size rather than one of sixteen, or the analysis is entirely within one projected CRS and you want the arithmetic to be plain Euclidean.
For a single city, either works. For a national pipeline that ingests data from multiple sources, H3's stable string key is usually worth the area variation.
Edge cases or notes
latlng_to_celltakes latitude first. So doescell_to_boundary's output. Shapely wants the opposite. Flip explicitly, both ways.- H3 v4 renamed everything.
geo_to_h3is nowlatlng_to_cell,h3_to_geo_boundaryiscell_to_boundary. Code from v3 will not run. - Use
cell_area, neveraverage_hexagon_area, to convert counts to densities. The difference was 12% in the example above. - A hexagon crossing the antimeridian produces a polygon spanning the globe. Split it, or work in a projected CRS.
sjoinwithpredicate="within"drops points exactly on a boundary in rare floating-point cases. Check the binned total against the input total.- Empty cells are absent, not zero. If the analysis needs zeros — most neighbourhood statistics do — reindex against all cells in the extent.
- Circumradius, not width. In a projected hex grid,
sizeis centre-to-vertex; the flat-to-flat width issize * √3. - 12 pentagons exist per resolution. They have five neighbours, not six, and break any code that assumes six.
Internal links
- The modifiable areal unit problem explained — what hexagons do and do not fix
- Kernel density explained — the continuous alternative to binning
- How to count points in polygons with GeoPandas — the underlying join
- How to find hotspots with Getis-Ord Gi* in Python — what the hex grid feeds into
- Spatial weights warn about islands — neighbours on a hex grid
- How to aggregate spatial data by region in GeoPandas — aggregating to irregular units instead
- Points plot in the ocean off Africa — the coordinate-order failure H3 invites
- How to calculate area and distance in GeoPandas correctly — why cell area needs a projection
FAQ
Why hexagons instead of squares?
Every hexagon has six neighbours at the same distance; a square has four at d and four at d√2. That difference biases any neighbourhood operation along the grid axes.
Which H3 resolution should I use?
Pick from the area ladder: resolution 9 is about 0.1 km², resolution 10 about 0.015 km². The steps are roughly seven-fold, so there is no fine control — use a projected grid if you need a specific size.
Are H3 cells equal-area?
No. At one resolution the area varies about 20% globally, though only about 0.4% within a single city. Always use cell_area on your own cells rather than the published average.
Why are my hexagons in the wrong place?
Coordinate order. latlng_to_cell takes latitude first, and cell_to_boundary returns (lat, lng) while Shapely expects (x, y). Flip in both directions.
Do hexagons solve the modifiable areal unit problem?
No. They remove the directional bias of a square grid. Size and origin are still arbitrary, and the results still move when you change them.
How do I change scale without re-binning?
H3 is hierarchical: h3.cell_to_parent(cell, coarser_resolution) is a string operation, no geometry required.
Should I use H3 or a projected grid?
H3 for a stable join key, multi-region work, or cheap scale changes. A projected grid when cell area must be exactly equal or you need a specific cell size.