How to Use H3 in DuckDB for Grid Aggregation at Scale
Problem statement
You have millions of points and want hexagon counts: per resolution, with areas and densities, as polygons a map can draw. Doing it in Python means a loop over h3.latlng_to_cell, a pandas groupby, a second loop for areas, a third for boundaries and a GeoDataFrame at the end.
DuckDB's community h3 extension moves all of that into one SQL statement. Measured on 13,464,117 GeoNames points read straight from Parquet:
resolution occupied cells median per cell time
3 17,210 124 0.47 s
5 451,015 9 0.60 s
7 5,798,146 1 0.89 s
9 12,159,083 1 1.30 s
The surprising measurement is not the speed. It is that the obvious shortcut — index once at a fine resolution, then roll up with h3_cell_to_parent — is ten times faster than re-indexing and puts 6.51% of the points in a different cell.
Quick answer
install h3 from community;
load h3;
select h3_latlng_to_cell(lat, lon, 7) as cell, -- latitude first
count(*) as n
from read_parquet('geonames.parquet')
where lat is not null
group by 1
order by n desc;On the full 13.46 million rows that returned 5,798,146 cells in 0.89 s. The cell is a UBIGINT; convert to the familiar string form only for display, with h3_h3_to_string(cell).
Step-by-step solution
1. Install and load the extension
import duckdb
con = duckdb.connect()
con.execute("install h3 from community; load h3;")
print(con.execute("select extension_version from duckdb_extensions() "
"where extension_name = 'h3'").fetchone())With DuckDB 1.5.5 this reported v1.5.5 and registered 76 functions whose names start with h3_. install downloads once into ~/.duckdb/extensions/<version>/; load runs per connection. Community extensions are built by the community repository rather than the core team, so from community is required.
2. Get the argument order right
h3_latlng_to_cell(lat, lng, res) — latitude first, like h3-py, and the opposite of st_point(x, y). Reversing the arguments does not raise: measured, it returned 897b8172ed3ffff for London instead of 89194ad14c3ffff, a cell 7,491 km away. If the input is a geometry column, extract explicitly:
select h3_latlng_to_cell(st_y(geom), st_x(geom), 7) from places;The result matched h3-py exactly: 617438095025111039 as an integer, 89194ad14c3ffff as a string, for the same point in both libraries.
3. Choose the resolution from the distribution, not the map
The table above is the decision. At resolution 5 the median occupied cell held 9 points; at 7 it held 1. A hexagon map where most cells hold one point shows where records exist, not how dense they are. Run the aggregation at two or three resolutions and read the median before drawing anything.
4. Keep the cell as an integer
The string form is for people. Measured on the same 13.46 million resolution-9 cells, written to Parquet with zstd and then grouped:
stored as file size group by
UBIGINT 41.4 MB 0.20 s
VARCHAR 57.7 MB 0.65 s
Integers are 28% smaller on disk and 3.3× faster to group. Grouping by h3_latlng_to_cell_string directly was slower too: 1.16 s against 0.89 s at resolution 7.
5. Add area and density in the same query
select cell, n,
h3_cell_area(cell, 'km^2') as km2,
n / h3_cell_area(cell, 'km^2') as per_km2
from (select h3_latlng_to_cell(lat, lon, 5) as cell, count(*) as n
from read_parquet('geonames.parquet') group by 1);The unit argument is required. Without it:
Binder Error: No function matches the given name and argument types 'h3_cell_area(UBIGINT)'.
The busiest resolution-5 cell held 6,716 points in 172.5 km² — 38.9 per km² — while the next held 4,718 in 245.8 km². Resolution-5 cells range from 127.8 to 305.1 km², so a count map and a density map rank the same cells differently.
6. Roll up only when an approximate answer is acceptable
Once a resolution-9 column is materialised, deriving resolution 5 from it is very fast:
h3_cell_to_parent(h9, 5), grouped 0.06 s 451,684 cells
h3_latlng_to_cell(lat, lon, 5), grouped 0.53 s 451,015 cells
The two answers differ. Measured over every row, 876,010 of 13,464,117 points (6.51%) have a resolution-9 parent that is not the resolution-5 cell containing the point, and 250,742 resolution-5 cells ended up with a different count. H3 children do not nest exactly inside their parents; Example 2 measures this on your own data.
7. Export polygons as GeoParquet
install spatial; load spatial;
copy (
select h3_h3_to_string(cell) as h3, n,
st_geomfromtext(h3_cell_to_boundary_wkt(cell)) as geom
from (select h3_latlng_to_cell(lat, lon, 5) as cell, count(*) as n
from read_parquet('geonames.parquet') group by 1)
) to 'h3_res5.parquet' (format parquet);451,015 hexagons in 1.41 s, 49.0 MB, read back by GeoPandas with CRS OGC:CRS84. Build the polygons after aggregating: 451,015 boundaries are cheap; 13 million are not.
8. Check the export for hexagons that wrap the world
57 of those 451,015 polygons were wider than 180° of longitude — cells straddling the antimeridian, drawn the long way round. They are a tiny fraction of the rows and would dominate any map of the Pacific. Example 3 flags them in the same statement that writes the file.
Code examples
Example 1 — counts, area and density in one call
import time
import duckdb
def h3_connection(threads=6):
con = duckdb.connect()
con.execute(f"set threads = {threads}")
con.execute("set enable_progress_bar = false")
con.execute("install h3 from community; load h3;")
con.execute("install spatial; load spatial;")
return con
def h3_counts(con, source, res, lat="lat", lng="lon", where="true"):
"""Points per H3 cell, with area and density, computed entirely in SQL."""
started = time.perf_counter()
frame = con.execute(f"""
select h3_h3_to_string(cell) as h3,
n,
h3_cell_area(cell, 'km^2') as km2,
n / h3_cell_area(cell, 'km^2') as per_km2
from (
select h3_latlng_to_cell({lat}, {lng}, {res}) as cell, count(*) as n
from {source}
where {lat} is not null and {lng} is not null and {where}
group by 1
)
""").df()
print(f"res {res}: {len(frame):,} cells in {time.perf_counter() - started:.2f}s, "
f"median {frame['n'].median():.0f} points per occupied cell")
return frame>>> con = h3_connection()
>>> f = h3_counts(con, "read_parquet('geonames.parquet')", 5)
res 5: 451,015 cells in 0.82s, median 9 points per occupied cell
>>> f7 = h3_counts(con, "read_parquet('geonames.parquet')", 7, where="feature_class = 'P'")
res 7: 2,985,431 cells in 2.23s, median 1 points per occupied cell
The timing includes building the pandas frame, which is why it is slower than the bare aggregate. Keep the result in DuckDB, or write it to a file, if the frame would be large.
Example 2 — measure what a roll-up would change
def rollup_disagreement(con, source, fine, coarse, lat="lat", lng="lon"):
"""How many points change coarse cell if you roll up instead of re-indexing."""
moved, total = con.execute(f"""
select count(*) filter (
where h3_cell_to_parent(h3_latlng_to_cell({lat}, {lng}, {fine}), {coarse})
<> h3_latlng_to_cell({lat}, {lng}, {coarse})),
count(*)
from {source}
where {lat} is not null
""").fetchone()
print(f"res {fine} -> {coarse}: {moved:,} of {total:,} points ({100 * moved / total:.2f}%) "
f"belong to a different res-{coarse} cell than their res-{fine} parent")
return moved / totalres 9 -> 5: 876,010 of 13,464,117 points (6.51%) belong to a different res-5 cell than their res-9 parent
res 8 -> 7: 961,507 of 13,464,117 points (7.14%) belong to a different res-7 cell than their res-8 parent
res 6 -> 5: 963,235 of 13,464,117 points (7.15%) belong to a different res-5 cell than their res-6 parent
One level of roll-up moves about 7% of points. Four levels at once moved slightly fewer, 6.51% — the errors at successive levels do not simply add up. If downstream work joins to data indexed with h3_latlng_to_cell at the coarse resolution, roll-ups will not line up with it.
Example 3 — export hexagons and flag the antimeridian
import os
import geopandas as gpd
def export_hexagons(con, counts_sql, path):
"""Write hexagon polygons to GeoParquet and flag the ones that wrap the map."""
con.execute(f"""
copy (
select h3, n, geom,
st_xmax(geom) - st_xmin(geom) > 180 as crosses_180
from (
select h3_h3_to_string(cell) as h3, n,
st_geomfromtext(h3_cell_to_boundary_wkt(cell)) as geom
from ({counts_sql})
)
) to '{path}' (format parquet)
""")
frame = gpd.read_parquet(path)
print(f"{len(frame):,} hexagons, {int(frame['crosses_180'].sum())} cross the antimeridian, "
f"{os.path.getsize(path) / 1e6:.1f} MB, crs {frame.crs}")
return frame451,015 hexagons, 57 cross the antimeridian, 49.0 MB, crs OGC:CRS84
counts_sql must return columns named cell and n. The flag lets a map layer filter or repair those rows rather than discovering them as a band across the Pacific.
Explanation
Why the whole pipeline belongs in SQL
The index computation itself is not slow in Python: h3-py indexed 1,000,000 points in 0.75 s with string cells and 0.49 s with the numpy_int API. The cost is everything around it — building a list of strings, hashing them in pandas, computing areas in another loop, building shapely polygons.
In DuckDB the cell is computed, hashed and grouped inside one vectorised plan across six threads, and only the aggregated rows reach Python. The whole measured script, which also materialised a 13.46-million-row table, peaked at 2.48 GB of resident memory.
Why integer cells are faster
An H3 index is a 64-bit integer; the hexadecimal string is a 15-character rendering of it. Grouping by the integer hashes 8 bytes and compares one machine word. Grouping by the string hashes and compares variable-length data, and Parquet stores it less compactly. The measured 3.3× difference is that overhead, repeated 13.46 million times.
Why rolling up is not the same as re-indexing
H3 uses an aperture of 7: each cell has seven children, arranged so that the children's union approximates the parent. Approximately. The outer six children each stick out of the parent on one side and leave a gap on another, so a point near the edge of a coarse cell can sit in a child whose parent is the neighbouring cell.
The fraction is stable: 7.14–7.15% for a single level, measured twice at different resolutions. It matches a separate geometric measurement: across 500 random resolution-5 cells, 7.14% of the seven children's combined area lay outside the parent. Points are spread evenly enough that the share of misplaced points equals the share of misplaced area. Rolling up is the right operation when you want a strict hierarchy; re-indexing is right when you want each point in the coarse cell that contains it.
Why H3 costs more than a rectangular grid
A floor(lon / 0.1), floor(lat / 0.1) grid over the same rows took 0.14 s. H3 has to convert each coordinate to a 3D vector, find the nearest icosahedron face, project onto that face and then walk the resolution digits. That is several times the arithmetic of a floor, and it buys equal-ish areas, uniform neighbours and a global index — the reasons to choose hexagons at all.
Edge cases or notes
install h3 from communityneeds network access the first time. Offline machines and containers need the extension file baked in, like the spatial extension.- Every DuckDB upgrade needs a reinstall, because extensions are cached per DuckDB version.
h3_cell_arearequires a unit —'km^2','m^2'or'rads^2'. Measured,'m^2'returned exactly a million times'km^2'.- Null coordinates should be filtered before indexing rather than left to form a group.
_stringvariants exist for most functions (h3_latlng_to_cell_string,h3_polygon_wkt_to_cells_string); use them at the edges only.- Joining to cells from h3-py works on either form:
h3_string_to_h3converts strings to the integers DuckDB uses. - Polygon fills are available in SQL through
h3_polygon_wkt_to_cells, with the same centre-containment rules as h3-py. - Very high resolutions make almost every point its own group; resolution 9 produced 12.16 million groups from 13.46 million rows, a hash table nearly the size of the input.
Internal links
- How to aggregate millions of points into a grid with DuckDB — the rectangular-grid version this builds on
- How to install and load the DuckDB spatial extension — the same install story, and what breaks
- Choosing an H3 resolution: cell size, counts and what each level can show — reading the median column
- H3 hierarchy explained: parents, children and why they do not nest exactly — the mechanism behind the 6.51%
- How to aggregate points into H3 hexagons and map them — the Python route for smaller data
- Fixing H3 cells in the wrong place: latitude and longitude order — the silent argument swap
- Fixing H3 hexagons that stretch across the map at the antimeridian — what to do with the 57 flagged rows
- DuckDB's memory model: streaming, spilling and limits — why fine resolutions cost memory
FAQ
How do I install the H3 extension in DuckDB?
Run install h3 from community, then load h3 in each connection. The install downloads once per DuckDB version; the load is per session.
Is DuckDB faster than h3-py for aggregation?
For the whole job, yes, because indexing, grouping, areas and boundaries stay in one engine. The index call alone is not the bottleneck: h3-py indexed a million points in 0.75 s, and DuckDB aggregated all 13.46 million at resolution 7 in 0.89 s.
Should I store H3 cells as strings or integers?
Integers. On 13.46 million resolution-9 cells the UBIGINT Parquet file was 41.4 MB against 57.7 MB, and grouping took 0.20 s against 0.65 s.
Can I compute a fine resolution once and roll it up?
Only if an approximate hierarchy is acceptable. Rolling resolution 9 up to 5 placed 6.51% of points in a different cell from indexing them at resolution 5 directly.
Why does h3_cell_area fail with a binder error?
It needs a unit argument. Call h3_cell_area(cell, 'km^2'); the single-argument form does not exist.
How do I get hexagon polygons out of DuckDB?
Use h3_cell_to_boundary_wkt with st_geomfromtext and COPY to Parquet. Check the width of each polygon, because cells on the antimeridian come out spanning the whole map.