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,785 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,713 of 30,785 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))

30785

235,876 points became 30,785 hexagons in 0.43 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.

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,876 places with population of 500 or more) and the full 13,472,195-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 9 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  30785.0     30785.0     30785.0
mean       8.0    149403.0     97273.0
std       17.0    805176.0    571422.0
min        1.0         0.0         0.0
25%        1.0      3338.0      2446.0
50%        2.0     19425.0     12108.0
75%        6.0     77486.0     44039.0
max      396.0  47982043.0  24874500.0

The group-by itself took 0.010 s. A median of 19,425 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 987 cells holding 50,000–60,000 people have densities from 23.5 to 55.0 people per km², and their order by density correlates with their order by population at only 0.453.

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,785 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
import numpy as np

np.random.seed(0)                      # FisherJenksSampled classifies a random 30% sample
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 [30713, 49, 15, 6, 0, 1, 1]
quantiles              counts [4400, 4397, 4397, 4397, 4398, 4398, 4398]
fisher-jenks sample    counts [25257, 3579, 1121, 442, 168, 121, 97]

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")
print(f"{int(hexes['n'].sum()):,} points -> {len(hexes):,} cells")

13,472,195 points -> 451,134 cells

All 13,472,195 GeoNames points aggregated to 451,134 cells in 0.77 s for the bare group-by, with four threads. The whole of step 8 — areas, boundaries and the GeoDataFrame — took 4.2 s (Example 4). 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°.

Code examples

The data used on this page

GeoNames' cities500 (every populated place of 500 people or more) and the full allCountries dump, both downloaded from download.geonames.org on 28 September 2026 and written to Parquet with lat/lon columns. GeoNames is edited daily, so a fresh download will differ from these counts in the last digits. The France outline is Natural Earth's 1:10m Admin 0 countries.

import os
import zipfile

import pandas as pd
import requests

GEONAMES = ["geonameid", "name", "asciiname", "alternatenames", "lat", "lon", "feature_class",
            "feature_code", "country_code", "cc2", "admin1", "admin2", "admin3", "admin4",
            "population", "elevation", "dem", "timezone", "modified"]

def fetch(url, path):
    if not os.path.exists(path):
        with requests.get(url, stream=True, timeout=600) as reply:
            reply.raise_for_status()
            with open(path, "wb") as out:
                for chunk in reply.iter_content(1 << 20):
                    out.write(chunk)

if not os.path.exists("cities500.parquet"):
    fetch("https://download.geonames.org/export/dump/cities500.zip", "cities500.zip")
    with zipfile.ZipFile("cities500.zip") as z, z.open("cities500.txt") as f:
        table = pd.read_csv(f, sep="\t", names=GEONAMES, quoting=3, keep_default_na=False,
                            dtype={"country_code": str, "cc2": str, "admin1": str, "admin2": str,
                                   "admin3": str, "admin4": str})
    table.to_parquet("cities500.parquet")

if not os.path.exists("geonames_all.parquet"):
    fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")
    with zipfile.ZipFile("allCountries.zip") as z, z.open("allCountries.txt") as f:
        table = pd.read_csv(f, sep="\t", names=GEONAMES, quoting=3, keep_default_na=False,
                            usecols=["geonameid", "lat", "lon", "feature_class"], dtype={"feature_class": str})
    table.to_parquet("geonames_all.parquet")

fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries.zip",
      "ne_10m_admin_0_countries.zip")
print(len(pd.read_parquet("cities500.parquet")), len(pd.read_parquet("geonames_all.parquet", columns=["lat"])))

235876 13472195

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,419 cells in 0.23 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.32 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,419 resolution-5 hexagons took 3.1 s at 150 dpi without outlines and 3.4 s with white ones, the fastest of three runs each. 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].

Example 4 — the timings on this page

import os
import time

import duckdb
import matplotlib
matplotlib.use("Agg")


def fastest(f, runs=3):
    best = None
    for _ in range(runs):
        t = time.perf_counter()
        result = f()
        elapsed = time.perf_counter() - t
        best = elapsed if best is None or elapsed < best else best
    return best, result


print(f"load average {os.getloadavg()[0]:.1f} on {os.cpu_count()} threads")


def quick_answer():
    p = pd.read_parquet("cities500.parquet")
    p["cell"] = [h3.latlng_to_cell(a, b, 4) for a, b in zip(p["lat"].tolist(), p["lon"].tolist())]
    h = p.groupby("cell").agg(n=("name", "size"), population=("population", "sum"))
    return gpd.GeoDataFrame(h, geometry=[Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(c)])
                                         for c in h.index], crs="EPSG:4326")

t, h = fastest(quick_answer)
print(f"quick answer: {len(places):,} points -> {len(h):,} cells in {t:.2f} s")

t_list, _ = fastest(lambda: [h3.latlng_to_cell(a, b, 4)
                             for a, b in zip(places["lat"].tolist(), places["lon"].tolist())])
t_apply, _ = fastest(lambda: places.apply(lambda r: h3.latlng_to_cell(r["lat"], r["lon"], 4), axis=1), runs=1)
print(f"assignment: comprehension {t_list:.2f} s, DataFrame.apply {t_apply:.2f} s ({t_apply / t_list:.0f}x)")

t, _ = fastest(lambda: places.groupby("cell").agg(n=("geonameid", "size"), population=("population", "sum"),
                                                  largest=("population", "max")))
print(f"group-by: {t:.3f} s")

t, a5 = fastest(lambda: aggregate_to_h3(places, 5, sums=["population"], maxes=["population"]))
print(f"Example 1: {len(a5):,} cells in {t:.2f} s")

fr6 = aggregate_to_h3(fr, 6)
t, fr_gdf = fastest(lambda: hexes_to_gdf(fr6, region=mainland, res=6))
print(f"Example 2: {len(fr_gdf):,} polygons, {int((fr_gdf['n'] == 0).sum()):,} empty, in {t:.2f} s")

world = hexes_to_gdf(a5)
for edge in (None, "white"):
    t, _ = fastest(lambda: hex_map(world, "population_per_km2", "world.png", edge=edge))
    print(f"Example 3, edge={edge}: {t:.1f} s, {os.path.getsize('world.png') / 1e6:.2f} MB")
print("France classes:", hex_map(fr_gdf, "n", "france.png"))

con = duckdb.connect()
con.execute("load h3; set threads = 4")
t, n = fastest(lambda: con.execute("""select count(*) from (select h3_latlng_to_cell(lat, lon, 5), count(*)
                                      from read_parquet('geonames_all.parquet') group by 1)""").fetchone()[0])
print(f"DuckDB, 4 threads: {n:,} cells, bare group-by {t:.2f} s")


def step_8():
    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()
    return gpd.GeoDataFrame(frame.drop(columns="wkt"), geometry=gpd.GeoSeries.from_wkt(frame["wkt"]),
                            crs="EPSG:4326")

t, _ = fastest(step_8)
print(f"DuckDB, 4 threads: step 8 in full, with areas, boundaries and the GeoDataFrame, {t:.2f} s")

load average 1.9 on 24 threads
quick answer: 235,876 points -> 30,785 cells in 0.43 s
assignment: comprehension 0.09 s, DataFrame.apply 0.83 s (9x)
group-by: 0.010 s
Example 1: 83,419 cells in 0.23 s
Example 2: 16,149 polygons, 7,295 empty, in 0.32 s
Example 3, edge=None: 3.1 s, 0.26 MB
Example 3, edge=white: 3.4 s, 0.22 MB
France classes: [7295, 5122, 2172, 1560]
DuckDB, 4 threads: 451,134 cells, bare group-by 0.77 s
DuckDB, 4 threads: step 8 in full, with areas, boundaries and the GeoDataFrame, 4.22 s

Example 5 — the other numbers on this page

from pyproj import Geod
from scipy.stats import spearmanr

print(f"cell area ratio at res 4: {agg['km2'].max() / agg['km2'].min():.2f}")
top = set(agg.nlargest(100, "population")["cell"]) & set(agg.nlargest(100, "density")["cell"])
print(f"Spearman population vs density {spearmanr(agg['population'], agg['density']).statistic:.4f}; "
      f"top-100 overlap {len(top)}")
band = agg[agg["population"].between(50_000, 60_000)]
print(f"50-60k band: {len(band)} cells, density {band['density'].min():.1f}-{band['density'].max():.1f} "
      f"({band['density'].max() / band['density'].min():.2f}x), "
      f"Spearman {spearmanr(band['population'], band['density']).statistic:.3f}")

offshore = counts[~counts.index.isin(land)]
print(f"offshore occupied cells in France: {len(offshore)}, holding {int(offshore.sum())} places")

geod = Geod(ellps="WGS84")
for c in agg["cell"]:
    if Polygon([(lng, lat) for lat, lng in h3.cell_to_boundary(c)]).bounds[2] < 0:
        continue
    ring = [(lng, lat) for lat, lng in h3.cell_to_boundary(c)]
    if max(x for x, _ in ring) - min(x for x, _ in ring) > 180:
        split = cell_polygon(c)
        whole = Polygon([(lng + 360 if lng < 0 else lng, lat) for lng, lat in ring])
        parts = sum(abs(geod.geometry_area_perimeter(g)[0]) for g in split.geoms) / 1e6
        print(f"antimeridian: {c}, parts {[tuple(round(v, 3) for v in g.bounds) for g in split.geoms]}, "
              f"{parts:.2f} km2 split, {abs(geod.geometry_area_perimeter(whole)[0]) / 1e6:.2f} km2 whole")

ei = mapclassify.EqualInterval(agg["population"], k=7)
print(f"equal-interval first break {ei.bins[0]:,.0f}; max {agg['population'].max():,}")

shanghai = h3.latlng_to_cell(31.22222, 121.45806, 5)
in_cell = places[[h3.latlng_to_cell(a, b, 5) == shanghai
                  for a, b in zip(places["lat"].tolist(), places["lon"].tolist())]]
print(f"{shanghai}: {in_cell['population'].sum():,} people, "
      f"{a5.loc[shanghai, 'population_per_km2']:,.0f} per km2; densest cell {a5['population_per_km2'].idxmax()}")
print(in_cell.nlargest(3, "population")[["name", "feature_code", "population"]].to_string(index=False))
pplx = places[places["feature_code"] == "PPLX"]
print(f"PPLX rows: {len(pplx):,}, {pplx['population'].sum() / 1e6:.0f} million people")

cells7 = pd.Series([h3.latlng_to_cell(a, b, 7) for a, b in zip(places["lat"].tolist(), places["lon"].tolist())])
cells7 = cells7.drop_duplicates().sample(20_000, random_state=0)
g7 = gpd.GeoDataFrame({"cell": cells7.values}, crs="EPSG:4326",
                      geometry=[cell_polygon(c) for c in cells7])
g7["lat"] = [abs(h3.cell_to_latlng(c)[0]) for c in g7["cell"]]
g7["ratio"] = g7.to_crs(3857).area / [h3.cell_area(c, "m^2") for c in g7["cell"]]
print(f"Web Mercator area / true area: equator (<5°) {g7.loc[g7['lat'] < 5, 'ratio'].mean():.2f}, "
      f"50-60° {g7.loc[g7['lat'].between(50, 60), 'ratio'].mean():.2f}, "
      f"max {g7['ratio'].max():.2f} at {g7.loc[g7['ratio'].idxmax(), 'lat']:.1f}°")

sample = hexes.sample(1000, random_state=0)
print("DuckDB area matches h3 on 1,000 cells:",
      all(abs(r.km2 - h3.cell_area(r.cell, "km^2")) < 1e-9 for r in sample.itertuples()))
wide = hexes.geometry.bounds.eval("maxx - minx") > 180
print(f"DuckDB polygons wider than 180°: {int(wide.sum())} of {len(hexes):,}")

cell area ratio at res 4: 2.38
Spearman population vs density 0.9982; top-100 overlap 91
50-60k band: 987 cells, density 23.5-55.0 (2.33x), Spearman 0.453
offshore occupied cells in France: 235, holding 335 places
antimeridian: 8417691ffffffff, parts [(178.965, 62.894, 180.0, 63.314), (-180.0, 63.068, -179.956, 63.137)], 1892.63 km2 split, 1892.62 km2 whole
equal-interval first break 6,854,578; max 47,982,043
8530995bfffffff: 45,458,645 people, 235,932 per km2; densest cell 8530995bfffffff
    name feature_code  population
Shanghai         PPLA    24874500
    Puxi         PPLX     6683712
  Pudong        PPLA4     5681512
PPLX rows: 10,546, 219 million people
Web Mercator area / true area: equator (<5°) 1.00, 50-60° 2.77, max 12.44 at 73.5°
DuckDB area matches h3 on 1,000 cells: True
DuckDB polygons wider than 180°: 57 of 451,134

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.453. 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,713 cells. Quantiles divide by rank instead, so every class holds about 4,398 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.

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,546 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.00 times the true area within 5° of the equator, 2.77 times at 50–60° and 12.44 times at 73.5°. 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,134 resolution-5 polygons came back wider than 180°.

Internal links

• How to assign points to H3 cells in Python — the assignment step, fast and with bad rows handled
• Choosing an H3 resolution — how to pick the level before aggregating
• How to turn H3 cells into a GeoDataFrame of polygons — faster polygon construction for large results
• Fixing an H3 map that misleads — unequal cell areas and empty cells on the finished map
• Fixing H3 hexagons that stretch across the map at the antimeridian — the stripe, in depth
• How to use H3 in DuckDB for grid aggregation at scale — the SQL route for very large inputs
• How to fill a polygon with H3 cells — the region fill behind the empty cells
• Choropleth classification explained — quantiles, equal intervals and natural breaks
• How to bin points into hexagons in Python — the planar alternative when a global grid is not needed
• How to build a species richness grid — counting distinct species per cell instead of points
• How to map forest loss hotspots with H3 — rasters to hexagons, normalised by baseline

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,876 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,713 of 30,785 cells in one class, and sampled natural breaks still put 25,257 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.47 million GeoNames points reduced to 451,134 cells in 0.77 s.