How to Use H3 in DuckDB for Grid Aggregation at Scale
On this page
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,472,195 GeoNames points read straight from Parquet (Example 4):
resolution occupied cells median per cell time
3 17,243 123 0.49 s
5 451,134 9 0.59 s
7 5,800,166 1 0.82 s
9 12,166,755 1 1.05 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.47 million rows that returned 5,800,166 cells in 0.82 s. The cell is a UBIGINT; convert to the familiar string form only for display, with h3_h3_to_string(cell).
What you get — from The data used on this page
13472195 points
Load the extension, mind the order, keep integers
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.47 million resolution-9 cells, written to Parquet with zstd and then grouped:
stored as file size group by
UBIGINT 41.4 MB 0.11 s
VARCHAR 57.7 MB 0.27 s
Integers are 28% smaller on disk and 2.5× faster to group. Grouping by h3_latlng_to_cell_string directly was slower too: 1.09 s against 0.82 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² (discrete global grids measures all of them), 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.04 s 451,801 cells
h3_latlng_to_cell(lat, lon, 5), grouped 0.45 s 451,134 cells
The two answers differ. Measured over every row, 876,538 of 13,472,195 points (6.51%) have a resolution-9 parent that is not the resolution-5 cell containing the point, and 250,777 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,134 hexagons, 48.9 MB, read back by GeoPandas with CRS OGC:CRS84 (Example 3). Build the polygons after aggregating: 451,134 boundaries are cheap; 13 million are not.
8. Check the export for hexagons that wrap the world
57 of those 451,134 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.
Counts, area and density in one SQL call
The data used on this page
The points are GeoNames' full allCountries dump — every feature, about 13.5 million, with its feature class — from download.geonames.org, rebuilt daily, so counts drift a little between runs. The setup writes latitude, longitude and feature class to geonames.parquet, which every query below reads directly.
import os
import zipfile
import duckdb
import pandas as pd
import requests
if not os.path.exists("allCountries.zip"):
with requests.get("https://download.geonames.org/export/dump/allCountries.zip",
stream=True, timeout=600) as reply:
reply.raise_for_status()
with open("allCountries.zip", "wb") as out:
for chunk in reply.iter_content(1 << 20):
out.write(chunk)
if not os.path.exists("geonames.parquet"):
with zipfile.ZipFile("allCountries.zip") as z, z.open("allCountries.txt") as f:
gn = pd.read_csv(f, sep="\t", header=None, usecols=[4, 5, 6], quoting=3,
names=["lat", "lon", "feature_class"], keep_default_na=False)
gn.to_parquet("geonames.parquet")
del gn
print(duckdb.sql("select count(*) from read_parquet('geonames.parquet')").fetchone()[0], "points")13472195 points
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 framecon = h3_connection()
f = h3_counts(con, "read_parquet('geonames.parquet')", 5)
f7 = h3_counts(con, "read_parquet('geonames.parquet')", 7, where="feature_class = 'P'")res 5: 451,134 cells in 0.84s, median 9 points per occupied cell
res 7: 2,986,958 cells in 2.13s, 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 / totalfor fine, coarse in ((9, 5), (8, 7), (6, 5)):
rollup_disagreement(con, "read_parquet('geonames.parquet')", fine, coarse)res 9 -> 5: 876,538 of 13,472,195 points (6.51%) belong to a different res-5 cell than their res-9 parent
res 8 -> 7: 962,112 of 13,472,195 points (7.14%) belong to a different res-7 cell than their res-8 parent
res 6 -> 5: 963,801 of 13,472,195 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 framehexes = export_hexagons(con, """select h3_latlng_to_cell(lat, lon, 5) as cell, count(*) as n
from read_parquet('geonames.parquet') group by 1""", "h3_res5.parquet")451,134 hexagons, 57 cross the antimeridian, 48.9 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.
Example 4 — the other measurements on this page
import h3
import numpy as np
def timed(sql, repeat=3):
best = None
for _ in range(repeat):
t = time.perf_counter()
out = con.execute(sql).fetchall()
best = min(best or 1e9, time.perf_counter() - t)
return out, best
src = "read_parquet('geonames.parquet')"
print(con.execute("select extension_version from duckdb_extensions() where extension_name = 'h3'").fetchone(),
duckdb.__version__,
con.execute("select count(distinct function_name) from duckdb_functions() where function_name like 'h3_%'").fetchone()[0],
"h3_ functions")
for res in (3, 5, 7, 9):
(row,), secs = timed(f"""select count(*), median(n) from
(select h3_latlng_to_cell(lat, lon, {res}) c, count(*) n from {src} where lat is not null group by 1)""")
print(f"res {res}: {row[0]:,} occupied cells, median {row[1]:.0f}, {secs:.2f} s")
right, wrong = con.execute("select h3_latlng_to_cell(51.5007, -0.1246, 9), h3_latlng_to_cell(-0.1246, 51.5007, 9)").fetchone()
print(f"London {h3.int_to_str(right)} ({right}); arguments reversed {h3.int_to_str(wrong)}, "
f"{h3.great_circle_distance(h3.cell_to_latlng(h3.int_to_str(right)), h3.cell_to_latlng(h3.int_to_str(wrong)), 'km'):,.0f} km away; "
f"h3-py agrees: {h3.latlng_to_cell(51.5007, -0.1246, 9) == h3.int_to_str(right)}")
con.execute(f"create or replace table cells9 as select h3_latlng_to_cell(lat, lon, 9) h9, lat, lon from {src}")
con.execute("copy (select h9 from cells9) to 'h9_int.parquet' (format parquet, compression zstd)")
con.execute("copy (select h3_h3_to_string(h9) h9 from cells9) to 'h9_str.parquet' (format parquet, compression zstd)")
for label, path in (("UBIGINT", "h9_int.parquet"), ("VARCHAR", "h9_str.parquet")):
_, secs = timed(f"select count(*) from (select h9, count(*) from read_parquet('{path}') group by 1)")
print(f"{label:8} {os.path.getsize(path) / 1e6:6.1f} MB group by {secs:.2f} s")
_, secs = timed(f"select count(*) from (select h3_latlng_to_cell_string(lat, lon, 7), count(*) from {src} group by 1)")
print(f"grouping by h3_latlng_to_cell_string at res 7: {secs:.2f} s")
busy = con.execute(f"""select n, h3_cell_area(cell, 'km^2') km2 from
(select h3_latlng_to_cell(lat, lon, 5) cell, count(*) n from {src} group by 1) order by n desc limit 2""").fetchall()
for n, km2 in busy:
print(f"res-5 cell: {n:,} points in {km2:.1f} km² = {n / km2:.1f} per km²")
try:
con.execute("select h3_cell_area(h3_latlng_to_cell(51.5, -0.12, 5))")
except duckdb.Error as e:
print(str(e).splitlines()[0])
m2, km2 = con.execute("select h3_cell_area(h3_latlng_to_cell(51.5, -0.12, 5), 'm^2'), "
"h3_cell_area(h3_latlng_to_cell(51.5, -0.12, 5), 'km^2')").fetchone()
print(f"'m^2' ÷ 'km^2' = {m2 / km2:,.1f}")
(rolled,), t_roll = timed("select count(*) from (select h3_cell_to_parent(h9, 5), count(*) from cells9 group by 1)")
(direct,), t_dir = timed("select count(*) from (select h3_latlng_to_cell(lat, lon, 5), count(*) from cells9 group by 1)")
changed = con.execute("""select count(*) from
(select h3_cell_to_parent(h9, 5) c, count(*) n from cells9 group by 1) a
full join (select h3_latlng_to_cell(lat, lon, 5) c, count(*) n from cells9 group by 1) b using (c)
where a.n is distinct from b.n""").fetchone()[0]
print(f"roll-up {t_roll:.2f} s, {rolled[0]:,} cells; re-index {t_dir:.2f} s, {direct[0]:,} cells; "
f"{changed:,} res-5 cells with a different count")
million = con.execute(f"select lat, lon from {src} using sample 1000000 rows (reservoir, 0)").df()
t = time.perf_counter()
[h3.latlng_to_cell(a, b, 9) for a, b in zip(million.lat, million.lon)]
t_str = time.perf_counter() - t
import h3.api.numpy_int as h3n
t = time.perf_counter()
[h3n.latlng_to_cell(a, b, 9) for a, b in zip(million.lat, million.lon)]
t_int = time.perf_counter() - t
print(f"h3-py, 1,000,000 points: {t_str:.2f} s strings, {t_int:.2f} s numpy_int")
_, secs = timed(f"select count(*) from (select floor(lon / 0.1), floor(lat / 0.1), count(*) from {src} group by 1, 2)")
print(f"0.1° floor grid: {secs:.2f} s; res-9 groups {con.execute('select count(distinct h9) from cells9').fetchone()[0]:,}; "
f"load average {os.getloadavg()[0]:.1f}")('v1.5.5',) 1.5.5 76 h3_ functions
res 3: 17,243 occupied cells, median 123, 0.49 s
res 5: 451,134 occupied cells, median 9, 0.59 s
res 7: 5,800,166 occupied cells, median 1, 0.82 s
res 9: 12,166,755 occupied cells, median 1, 1.05 s
London 89194ad14c3ffff (617438095025111039); arguments reversed 897b8172ed3ffff, 7,491 km away; h3-py agrees: True
UBIGINT 41.4 MB group by 0.11 s
VARCHAR 57.7 MB group by 0.27 s
grouping by h3_latlng_to_cell_string at res 7: 1.09 s
res-5 cell: 6,716 points in 172.5 km² = 38.9 per km²
res-5 cell: 4,718 points in 245.8 km² = 19.2 per km²
Binder Error: No function matches the given name and argument types 'h3_cell_area(UBIGINT)'. You might need to add explicit type casts.
'm^2' ÷ 'km^2' = 1,000,000.0
roll-up 0.04 s, 451,801 cells; re-index 0.45 s, 451,134 cells; 250,777 res-5 cells with a different count
h3-py, 1,000,000 points: 0.60 s strings, 0.46 s numpy_int
0.1° floor grid: 0.10 s; res-9 groups 12,166,755; load average 3.5
Why the whole pipeline belongs in SQL
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.60 s with string cells and 0.46 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.
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 2.5× difference is that overhead, repeated 13.47 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.10 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.
Network installs, reinstalls and unit arguments
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.17 million groups from 13.47 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.60 s, and DuckDB aggregated all 13.47 million at resolution 7 in 0.82 s.
Should I store H3 cells as strings or integers?
Integers. On 13.47 million resolution-9 cells the UBIGINT Parquet file was 41.4 MB against 57.7 MB, and grouping took 0.11 s against 0.27 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.