⚙️ How-To H3 & global grids Spatial statistics Maps & visualisation H3 & Global Grids · step 5 of 20 Updated

How to Aggregate Points into H3 Hexagons and Map Them

Problem statement

You have points, such as places, incidents, sensor readings or customers, and you want a hexagon map: one value per cell, coloured by density. Each step is short, and several can mislead without raising an error:

  • Cell areas are not equal. The 30,775 occupied resolution-4 cells in GeoNames' cities500 range from 896.6 km² to 2,136.0 km², a ratio of 2.38. Raw counts in two cells therefore compare different amounts of ground.
  • Empty cells are missing, not zero. Aggregating points gives rows only where points exist. Over mainland France at resolution 6, 45.8% of land cells have no place, and averaging only occupied cells nearly doubles the mean (1.74 against 0.94).
  • The default classification fails on skewed data. Equal intervals put 30,703 of 30,775 cells into one colour class.
  • One cell crosses the antimeridian and, drawn naively, becomes a 360°-wide stripe across the whole map.
  • Summing a gazetteer's populations can double count. The resolution-5 cell containing Shanghai sums to 45.5 million, because GeoNames lists the city and its districts as separate places.

This guide covers the whole pipeline, points to cells to aggregates to polygons to a map, with each of those problems handled where it arises.

Quick answer 

import geopandas as gpd
import h3
import pandas as pd
from shapely.geometry import Polygon

places = pd.read_parquet("cities500.parquet")
places["cell"] = [h3.latlng_to_cell(lat, lng, 4)
                  for lat, lng in zip(places["lat"].tolist(), places["lon"].tolist())]
hexes = places.groupby("cell").agg(n=("name", "size"), population=("population", "sum"))
hexes = gpd.GeoDataFrame(
    hexes,
    geometry=[Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(c)]) for c in hexes.index],
    crs="EPSG:4326",
)
print(len(hexes))
30775

235,735 points became 30,775 hexagons in 0.48 s, including the Parquet read and the polygon build. Before you map the result, divide by h3.cell_area, add the empty cells, choose a classification that suits skewed data, and split the one antimeridian cell. Each is a step below.

Flow from points through H3 cell assignment, group-by, area normalisation and polygon building to a hexagon map.
Polygons come last: the group-by works on 64-bit identifiers, and geometry is built only for the cells that survive.

Step-by-step solution

1. Load the points and choose a resolution

Any table with latitude and longitude columns in degrees works. The examples use GeoNames' cities500 (235,735 places with population of 500 or more) and the full 13,464,117-row GeoNames dump for scale. Pick the resolution for the pattern you want to show, not for the precision of the points. For a world map of settlements, resolution 4 (edge 26 km) or 5 (edge 9.9 km) is appropriate. Choosing an H3 resolution shows how to check the choice against your data.

2. Assign every point to a cell

Use the list comprehension from the quick answer. DataFrame.apply does the same job about 12 times more slowly, and argument order matters: latitude comes first. How to assign points to H3 cells covers speed, integer storage and bad coordinates.

3. Group by cell 

agg = (places.groupby("cell")
             .agg(n=("geonameid", "size"),
                  population=("population", "sum"),
                  largest=("population", "max"))
             .reset_index())
print(agg.describe().round(0))
             n  population     largest
count  30775.0     30775.0     30775.0
mean       8.0    149420.0     97289.0
std       17.0    805240.0    571487.0
min        1.0         0.0         0.0
25%        1.0      3334.0      2446.0
50%        2.0     19429.0     12108.0
75%        6.0     77482.0     44020.0
max      396.0  47982043.0  24874500.0

The group-by itself took 0.014 s. A median of 19,429 people against a maximum of nearly 48 million is the skew that every later step has to handle.

4. Divide by each cell's own area 

agg["km2"] = [h3.cell_area(c, "km^2") for c in agg["cell"]]
agg["density"] = agg["population"] / agg["km2"]
print(agg["km2"].agg(["min", "median", "max"]).round(1))
min        896.6
median    1793.0
max       2136.0
Name: km2, dtype: float64

Across all cells, normalising barely changes the ranking: the Spearman correlation between population and density is 0.9982, and 91 of the top 100 cells are the same on both measures. Among cells that are otherwise similar, it changes a great deal. The 983 cells holding 50,000–60,000 people have densities from 23.6 to 55.0 people per km², and their order by density correlates with their order by population at only 0.445.

5. Add the empty cells of the study area

A group-by produces rows only where there are points. To show where there are none, fill the region with cells and reindex:

countries = gpd.read_file("ne_10m_admin_0_countries.zip")
france = countries.loc[countries["ADMIN"] == "France", "geometry"].iloc[0]
mainland = max(france.geoms, key=lambda part: part.area)

fr = places[places["country_code"] == "FR"]
counts = pd.Series([h3.latlng_to_cell(a, b, 6)
                    for a, b in zip(fr["lat"].tolist(), fr["lon"].tolist())]).value_counts()
land = pd.Index(sorted(h3.geo_to_cells(mainland, 6)), name="cell")
full = counts.reindex(land, fill_value=0)
print(f"{len(land):,} land cells, {int((full == 0).sum()):,} empty ({(full == 0).mean():.1%})")
print(f"mean per land cell {full.mean():.2f}; mean per occupied cell {counts.mean():.2f}")
15,914 land cells, 7,295 empty (45.8%)
mean per land cell 0.94; mean per occupied cell 1.74

Reindexing also drops occupied cells whose centre is offshore: 235 cells holding 335 places along the coast and on small islands. Check that number before deciding whether that loss is acceptable.

6. Build polygons, splitting the antimeridian cell

h3.cell_to_boundary returns (lat, lng) pairs, so swap them for Shapely. A cell that straddles 180° has vertices near +180 and −180, and a naive polygon joins them the long way round:

from shapely import MultiPolygon, affinity
from shapely.geometry import Polygon, box


def cell_polygon(cell):
    ring = [(lng, lat) for lat, lng in h3.cell_to_boundary(cell)]
    poly = Polygon(ring)
    if poly.bounds[2] - poly.bounds[0] <= 180:
        return poly
    shifted = Polygon([(lng + 360 if lng < 0 else lng, lat) for lng, lat in ring])
    east = shifted.intersection(box(0, -90, 180, 90))
    west = affinity.translate(shifted.intersection(box(180, -90, 540, 90)), xoff=-360)
    return MultiPolygon([east, west])

Of the 30,775 occupied cells, exactly one crosses: 8417691ffffffff in Chukotka. Split, it becomes two parts with bounds (178.965, 62.894, 180.0, 63.314) and (-180.0, 63.068, -179.956, 63.137), and the total area is unchanged. The dedicated fix is in Fixing H3 hexagons that stretch across the map.

7. Classify for skewed data, then map in an equal-area projection 

import mapclassify

for name, cls in [("equal interval", mapclassify.EqualInterval(agg["population"], k=7)),
                  ("quantiles", mapclassify.Quantiles(agg["population"], k=7)),
                  ("fisher-jenks sample", mapclassify.FisherJenksSampled(agg["population"], k=7, pct=0.3))]:
    print(f"{name:22} counts {cls.counts.tolist()}")
equal interval         counts [30703, 49, 15, 6, 0, 1, 1]
quantiles              counts [4399, 4394, 4397, 4396, 4396, 4396, 4397]
fisher-jenks sample    counts [28574, 1653, 347, 104, 50, 26, 21]

Equal intervals give a single-colour map with a few dots. Natural breaks are not much better on data this skewed. Quantiles use every colour, at the cost of breaks the reader has to read from the legend. Plot in an equal-area projection such as Equal Earth (EPSG:8857), so that hexagons at 60° do not look larger than those at the equator.

8. Move the aggregation into DuckDB when the input is large 

import duckdb
import geopandas as gpd

con = duckdb.connect()
con.execute("install h3 from community; load h3")
frame = con.execute("""
    with counts as (
        select h3_latlng_to_cell(lat, lon, 5) as cell, count(*) as n
        from read_parquet('geonames_all.parquet')
        group by 1
    )
    select h3_h3_to_string(cell) as cell, n,
           h3_cell_area(cell, 'km^2') as km2,
           h3_cell_to_boundary_wkt(cell) as wkt
    from counts
""").df()
hexes = gpd.GeoDataFrame(frame.drop(columns="wkt"),
                         geometry=gpd.GeoSeries.from_wkt(frame["wkt"]), crs="EPSG:4326")

All 13,464,117 GeoNames points aggregated to 451,015 cells in 1.27 s for the bare group-by, with four threads on a shared machine. Adding areas and boundaries and building the GeoDataFrame took 4.5 s. DuckDB's h3_cell_area matched Python's h3.cell_area exactly on a 1,000-cell sample. Its boundary WKT has the same antimeridian problem: 57 of the polygons were wider than 180°.

Grid comparing how many of 30,775 hexagons fall into each of seven classes under equal intervals and under quantiles.
A class holding 99.8% of the cells is a legend entry, not information.

Code examples

Example 1 — aggregate with counts, sums, maxima and densities 

import h3
import pandas as pd


def aggregate_to_h3(df, res, sums=(), maxes=(), lat="lat", lon="lon"):
    """Count points per H3 cell, sum and max chosen columns, add area and density."""
    cells = [h3.latlng_to_cell(a, b, res) for a, b in zip(df[lat].tolist(), df[lon].tolist())]
    grouped = df.assign(cell=cells).groupby("cell")
    out = grouped.size().rename("n").to_frame()
    for col in sums:
        out[f"{col}_sum"] = grouped[col].sum()
    for col in maxes:
        out[f"{col}_max"] = grouped[col].max()
    out["km2"] = [h3.cell_area(c, "km^2") for c in out.index]
    out["n_per_km2"] = out["n"] / out["km2"]
    for col in sums:
        out[f"{col}_per_km2"] = out[f"{col}_sum"] / out["km2"]
    return out

Called as aggregate_to_h3(places, 5, sums=["population"], maxes=["population"]), it returned 83,391 cells in 0.30 s for cities500. The densest cell by population is 8530995bfffffff (Shanghai) at 235,932 people per km², and that number is inflated. The edge cases below explain why.

Example 2 — polygons for an aggregate, with the empty cells of a region 

import geopandas as gpd
import h3
import pandas as pd


def hexes_to_gdf(frame, region=None, res=None, fill=0):
    """Polygons for an aggregated frame indexed by cell; optionally add the empty cells of a region."""
    if region is not None:
        region_cells = pd.Index(sorted(h3.geo_to_cells(region, res)), name=frame.index.name)
        frame = frame.reindex(frame.index.union(region_cells))
        count_cols = [c for c in frame.columns if c == "n" or c.endswith(("_sum", "_per_km2"))]
        frame[count_cols] = frame[count_cols].fillna(fill)
        frame["km2"] = [h3.cell_area(c, "km^2") for c in frame.index]
    geoms = [cell_polygon(c) for c in frame.index]
    return gpd.GeoDataFrame(frame, geometry=geoms, crs="EPSG:4326")

With the France aggregate at resolution 6 and the mainland polygon as region, it returned 16,149 polygons, 7,295 of them empty, in 0.52 s. The union keeps the 235 offshore cells, which is why there are more rows than land cells. cell_polygon is the antimeridian-safe function from step 6.

Example 3 — a quantile hexagon map in an equal-area projection 

import mapclassify
import matplotlib.pyplot as plt


def hex_map(gdf, column, path, k=7, crs="EPSG:8857", cmap="viridis", edge=None):
    """Quantile hexagon map in an equal-area projection; returns the class counts."""
    projected = gdf.to_crs(crs)
    classes = mapclassify.Quantiles(projected[column], k=k)
    fig, ax = plt.subplots(figsize=(10, 5.5))
    projected.plot(column=column, ax=ax, cmap=cmap, scheme="quantiles", k=k,
                   linewidth=0 if edge is None else 0.2, edgecolor=edge or "face",
                   legend=True, legend_kwds={"loc": "lower left", "fontsize": 7})
    ax.set_axis_off()
    fig.savefig(path, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return classes.counts.tolist()

Rendering all 83,391 resolution-5 hexagons took between 5.3 and 5.9 s per map at 150 dpi over three runs. White outlines made no consistent difference to the time. The PNG was 0.26 MB without outlines and 0.22 MB with them. For the France map, which is mostly zeros and ones, the function returned four classes, not seven: [7295, 5122, 2172, 1560].

Explanation

Why normalising matters most where it seems not to

Across the whole world, population per cell runs from zero to 48 million, while cell area varies by a factor of 2.4. Dividing one by the other hardly moves the ranking, which is why the overall Spearman correlation is 0.998 and why skipping the step can look harmless.

A map, though, is read locally: this cell against its neighbour, this region against that one. Among cells of similar population, area is a large share of the remaining variation. In the 50,000–60,000 band, densities differ by a factor of 2.33 and the two orderings correlate at 0.445. That is the comparison readers actually make, and area alone decides it. The mechanics, and what else can mislead on a finished map, are in Fixing an H3 map that misleads.

Why the empty cells change the answer

An aggregate grouped from points is a list of places where something happened. Every statistic computed on it, whether mean, median or a quantile break, describes occupied cells only. Over mainland France, occupied cells average 1.74 places each, and all land cells average 0.94. The first figure answers "how many places does a cell with a place have", which is rarely the question.

Filling the region also changes the classification. With 7,295 zeros, quantile breaks collapse: mapclassify warns Not enough unique values in array to form 7 classes. Setting k to 4. That warning is correct. Map zeros as their own class, or classify only the non-zero cells and draw the empty ones in a neutral colour.

Why equal intervals fail on counts

Equal intervals divide the range, and the range of a count is set by its maximum, here one cell with 48 million people. Seven equal slices of 0–47,982,043 make the first break 6,854,578, so every cell below that falls into the first class: 30,703 cells. Quantiles divide by rank instead, so every class holds about 4,396 cells. Background on the trade-offs is in Choropleth classification explained.

Why one hexagon can ruin a world map

H3 cells are defined on the sphere, but Shapely polygons live on a flat plane where −180 and +180 are 360 units apart. A cell whose vertices sit at 179.9° and −179.9° becomes a polygon whose edges run the long way round, a thin band across every longitude. At resolution 4 only one occupied cell does this. At resolution 5 in DuckDB, 57 do. On a filled polygon map, a single one covers the whole width of the map. Splitting at the antimeridian keeps the area and gives two valid parts that any projection can draw.

Two schematic hexagon maps: occupied cells only, and the same region with empty cells filled as zero.
The left map cannot show where nothing is, and any mean computed from it describes occupied ground only.

Edge cases or notes

  • Gazetteer populations overlap. The Shanghai cell sums the city (24,874,500) and districts such as Puxi (6,683,712) and Pudong (5,681,512). GeoNames has 10,542 PPLX "section of populated place" rows in cities500 holding 219 million people. Filter feature codes before summing.
  • Filling uses cell centres. geo_to_cells keeps cells whose centre is inside the polygon, so coastal cells can fall out. On France that was 235 occupied cells.
  • Web Mercator distorts hexagon areas visually. On 20,000 occupied resolution-7 cells, EPSG:3857 areas were 1.02 times the true area near the equator, 2.65 times at 50–60° and 9.34 times at 70.9°. Map in an equal-area CRS.
  • Quantile breaks can repeat. Many tied values, zeros especially, make mapclassify reduce k and warn. Treat zero as its own class.
  • Counts per cell are not densities of points. Divide by h3.cell_area for a density. The table value from average_hexagon_area is the same for every cell and corrects nothing.
  • DuckDB's boundary WKT is not antimeridian-safe either. 57 of 451,015 resolution-5 polygons came back wider than 180°.

FAQ

How do I count points per H3 hexagon in Python?

Assign each point with h3.latlng_to_cell(lat, lng, res) in a list comprehension, then group by the cell column. For 235,735 places at resolution 4 that took under half a second, including building the polygons.

Should I divide H3 counts by cell area?

Yes, whenever cells will be compared. Resolution-4 cells vary from 896.6 to 2,136.0 km². Among cells with similar counts, area alone decided the order: densities spanned a factor of 2.33.

How do I show hexagons with no points?

Fill the study area with h3.geo_to_cells and reindex the aggregate on those cells with a fill value of zero. Over mainland France at resolution 6, 45.8% of land cells were empty.

Why does my H3 map have a stripe across it?

A cell that crosses the antimeridian was drawn as a polygon with vertices near +180 and −180. Split such cells at 180° before plotting. One occupied resolution-4 cell needed it.

Which classification works for hexagon counts?

Quantiles, or a scheme built for skewed data. Equal intervals put 30,703 of 30,775 cells in one class, and sampled natural breaks still put 28,574 there.

How do I aggregate tens of millions of points?

Do the group-by in DuckDB with the h3 extension and build polygons only for the result. All 13.46 million GeoNames points reduced to 451,015 cells in 1.27 s.