How to Bin Points into Hexagons in Python

On this page

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
from shapely.geometry import Polygon

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")

1637 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 by tens of percent globally	exactly equal
works across CRSs	yes, globally	one projection only
hierarchy	parent/child built in	none

The square grid has two kinds of neighbour. That difference biases every neighbourhood operation along the axes.

A resolution from the cell size, ids, then polygons

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"))
leeds = zones.to_crs("EPSG:4326").union_all()          # the city boundary
print(f"{len(summary)} occupied of {len(h3.geo_to_cells(leeds, 9)):,} cells inside Leeds")

1637 occupied of 5,936 cells inside Leeds

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}")

92,260–92,880 m² (spread 0.7%), nominal average 105,333

Two things worth knowing. Within one city the cells are near-identical — a 0.7% spread across Leeds, so local comparisons are safe. But the nominal average is 13–14% 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.

The same resolution, four places, a 28% range — and not in latitude order. Locally near-constant, globally not.

H3 and projected hex grids for Leeds crimes

The data used on this page

Every example runs on real data: the 488 Lower Layer Super Output Areas (LSOAs) of Leeds from the ONS, their Census 2021 population from Nomis, and the street-level crimes West Yorkshire Police published for June 2026 on data.police.uk. The binning examples use the crime points, incidents.

import io
import time

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

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query")
features, offset = [], 0
while True:                                   # the service pages at 1,000 features
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD,LSOA21NM",
                                     "outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700).set_index("LSOA21CD").sort_index()

population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(zones.index)}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = population

west, south, east, north = zones.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                  data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if reply.status_code != 429:
                break
            time.sleep(2 ** attempt)
        reply.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "lon": float(c["location"]["longitude"]),
                    "lat": float(c["location"]["latitude"])} for c in reply.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes.lon, crimes.lat), crs=4326).to_crs(27700)
incidents = gpd.sjoin(incidents, zones[["geometry"]], predicate="within").rename(columns={"index_right": "lsoa", "LSOA21CD": "lsoa"})
zones["crimes"] = incidents["lsoa"].value_counts().reindex(zones.index, fill_value=0)
zones["rate"] = 1000 * zones["crimes"] / zones["population"]
print(f"{len(zones)} LSOAs, {zones['population'].sum():,} residents, {len(incidents):,} crimes; "
      f"rate median {zones['rate'].median():.1f}, max {zones['rate'].max():.1f} per 1,000")

488 LSOAs, 811,957 residents, 9,229 crimes; rate median 8.4, max 174.4 per 1,000

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)
print(hexes.nlargest(3, "n")[["h3", "n", "per_km2"]].round(0).to_string(index=False))

res 9: 1,637 occupied cells, 9,229 points, max 186 per cell (2,007/km²)
             h3   n  per_km2
891942c5d6fffff 186   2007.0
891942c5d63ffff 135   1457.0
891942c5d6bffff 135   1457.0

Computing per_km2 from each cell's own area rather than the nominal average is what makes those densities correct. The difference is 13–14% 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,644 cells, 1,547 occupied
cell area 103,923–103,923 m² (exactly equal: True)
points binned: 9,229 of 9,229

points binned: 9,229 of 9,229 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.

H3 buys a global key at the cost of slightly unequal cells; a projected grid is exact but local.

Example 3 — using the H3 hierarchy to change scale for free

def resolution_sweep(points, boundary, resolutions=(7, 8, 9, 10)):
    rows = []
    for r in resolutions:
        hexes = h3_bin(points, r)
        in_area = len(h3.h3shape_to_cells_experimental(      # every cell touching the area, occupied or not
            h3.geo_to_h3shape(boundary), r, contain="overlap"))
        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) / in_area:.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, leeds).to_string(index=False))

res 7: 125 occupied cells, 9,229 points, max 1124 per cell (248/km²)
res 8: 482 occupied cells, 9,229 points, max 695 per cell (1,071/km²)
res 9: 1,637 occupied cells, 9,229 points, max 186 per cell (2,007/km²)
res 10: 3,338 occupied cells, 9,229 points, max 80 per cell (6,043/km²)
 res  cells  median_n  max_n empty_share top_decile_share
   7    125        23   1124         21%              52%
   8    482         9    695         49%              45%
   9   1637         3    186         74%              41%
  10   3338         2     80         92%              38%

top_decile_share is the concentration measure: at resolution 7 the busiest 10% of occupied cells hold 52% of crimes; at resolution 10 they hold 38%, while the share of Leeds' cells with no crime at all rises from 21% to 92%. Same data — and the direction is not the one most people guess, because at resolution 10 the median occupied cell holds two crimes. 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 involved

Six equal neighbours, unequal areas and MAUP

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:

for place, lat, lng in (("equator", 0.0, 0.0), ("Leeds", 53.80, -1.55), ("Oslo", 59.91, 10.75),
                        ("Svalbard", 78.22, 15.65)):
    print(f"  res 9 at {place:<9} {h3.cell_area(h3.latlng_to_cell(lat, lng, 9), 'm^2'):>9,.0f} m²")

  res 9 at equator      78,386 m²
  res 9 at Leeds        92,676 m²
  res 9 at Oslo         75,425 m²
  res 9 at Svalbard     96,198 m²

A 28% range at one resolution, and not a smooth function of latitude: Oslo's cells are smaller than the equator's, Svalbard's the largest of the four. Within a single city the variation is negligible — 0.7% across Leeds — 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 moves by 14 points between resolution 7 and 10 on identical data, and the empty share by 71. 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.

Quote a concentration figure without its resolution and you have said nothing.

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.

Latitude first, v4 names and empty cells

latlng_to_cell takes latitude first. So does cell_to_boundary's output. Shapely wants the opposite. Flip explicitly, both ways.
H3 v4 renamed everything. geo_to_h3 is now latlng_to_cell, h3_to_geo_boundary is cell_to_boundary. Code from v3 will not run.
Use cell_area, never average_hexagon_area, to convert counts to densities. The difference was 13–14% in the example above.
A hexagon crossing the antimeridian produces a polygon spanning the globe. Split it, or work in a projected CRS.
sjoin with predicate="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, size is centre-to-vertex; the flat-to-flat width is $\text{size} \times \sqrt{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
How to aggregate points into H3 hexagons and map them — hexagons on a global grid instead of a local one
Choosing an H3 resolution — picking the cell size from the data

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 varied by 28% across four places measured here, though only 0.7% within Leeds. 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.