How to Join Two Point Datasets on an H3 Index Instead of a Spatial Join
Problem statement
You have two point layers and want every pair within some distance: hotels within 500 m of a railway station, sensors near a road junction, customers near a store. The textbook tool is a spatial join with a distance predicate, which needs geometries and a projected CRS whose metres are true enough. An H3 index looks like a shortcut. Give both sides a cell column and merge on it — an ordinary equality join that runs in pandas, SQL or Spark with no geometry at all.
Done naively, that shortcut is badly wrong. On 102,602 GeoNames hotels and 33,485 stations in Europe, a same-cell join at res 9 found 2,216 of the 13,931 true pairs within 500 m — 15.9%. Pairs that straddle a cell edge share no key, and at any useful resolution most do.
The fix keeps the merge and changes what is merged: expand one side into a disk of cells, join, and refine by true distance. With k chosen by a rule verified below, the same join returned 13,931 of 13,931 pairs, with no false ones.
Quick answer
import math
import numpy as np
import h3
def haversine_m(lat1, lon1, lat2, lon2):
p1, p2 = np.radians(lat1), np.radians(lat2)
a = np.sin((p2 - p1) / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(np.radians(lon2 - lon1) / 2) ** 2
return 2 * 6_371_008.8 * np.arcsin(np.sqrt(a))
res, distance = 9, 500
k = math.ceil(distance / (1.5 * h3.average_hexagon_edge_length(res, "m"))) + 1 # 3
a = hotels.assign(cell=[h3.latlng_to_cell(y, x, res) for y, x in zip(hotels.lat, hotels.lon)])
b = stations.assign(cell=[h3.grid_disk(h3.latlng_to_cell(y, x, res), k)
for y, x in zip(stations.lat, stations.lon)]).explode("cell")
pairs = a.merge(b, on="cell", suffixes=("_hotel", "_station"))
pairs = pairs[haversine_m(pairs.lat_hotel, pairs.lon_hotel, pairs.lat_station, pairs.lon_station) <= distance]
print(k, len(pairs))3 13931
The rule k = ceil(d / (1.5 × edge)) + 1 produced zero misses on 140,000 adversarial pairs placed at exactly the join distance, across seven resolution and distance combinations and around the pentagons. The refinement removes the candidates that were only close in cell terms.
Step-by-step solution
1. Load both sides and build a ground truth once
import duckdb
con = duckdb.connect()
con.execute("set threads = 4")
window = "lat between 34 and 72 and lon between -25 and 45"
hotels = con.execute(f"select geonameid as id, lat, lon from read_parquet('geonames_all.parquet') "
f"where feature_code = 'HTL' and {window}").df()
stations = con.execute(f"select geonameid as id, lat, lon from read_parquet('geonames_all.parquet') "
f"where feature_code in ('RSTN', 'RSTP') and {window}").df()That gives 102,602 hotels and 33,485 stations. An exact answer on the sphere is cheap at this size: put both sides on the unit sphere as x, y, z and ask a KD-tree for every pair within the chord that corresponds to 500 m.
from scipy.spatial import cKDTree
def unit_xyz(lat, lon):
la, lo = np.radians(lat), np.radians(lon)
return np.column_stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)])
chord = 2 * math.sin(distance / 6_371_008.8 / 2)
tree_h, tree_s = cKDTree(unit_xyz(hotels.lat, hotels.lon)), cKDTree(unit_xyz(stations.lat, stations.lon))
truth = {(hotels.id.iat[i], stations.id.iat[j])
for i, js in enumerate(tree_h.query_ball_tree(tree_s, chord)) for j in js}
print(f"{len(truth):,} pairs within {distance} m")13,931 pairs within 500 m
It took 0.33 s. Every other method below is scored against these 13,931 pairs.
2. See why a same-cell join is not enough
res 7: 50,163 candidates, 11,523 within 500 m, recall 82.7%
res 8: 11,336 candidates, 8,297 within 500 m, recall 59.6%
res 9: 2,216 candidates, 2,216 within 500 m, recall 15.9%
A finer cell misses more pairs because a 500 m separation crosses more edges. A coarser cell misses fewer but floods the join with candidates: at res 7, 38,640 of 50,163 candidate pairs were further apart than 500 m, and 17.3% of the true pairs were still missed. No resolution makes a same-cell join correct.
3. Choose k from the distance and the resolution
for r in (7, 8, 9, 10):
edge = h3.average_hexagon_edge_length(r, "m")
print(r, round(edge), [math.ceil(d / (1.5 * edge)) + 1 for d in (100, 250, 500, 1000, 5000)])7 1406 [2, 2, 2, 2, 4]
8 531 [2, 2, 2, 3, 8]
9 201 [2, 2, 3, 5, 18]
10 76 [2, 4, 6, 10, 45]
Each ring of a disk extends its guaranteed reach by about 1.4 edge lengths, so dividing by 1.5 rounds towards more rings; the extra ring covers the offset of each point from its own cell's centre. Tested with 20,000 random pairs placed at exactly the join distance for each of seven combinations, the rule's k missed nothing, while k − 1 missed up to 75 pairs (res 8 at 2 km). Around the twelve res-8 pentagons, 24,000 more pairs at 499 m produced no misses either.
Pick a resolution that keeps k at 2 or 3. Finer cells with a larger k give the same answer with far more keys: 5,000 m at res 10 needs k = 45, or 6,211 cells per point.
4. Expand the smaller side into disks and merge
grid_disk on each station, then explode so each cell becomes a row. At res 9 and k = 3 every disk holds 37 cells, so the 33,485 stations became 1,238,945 keys. Expanding the hotels instead would have made 3,796,274.
A disk never repeats a cell, so the merge — a plain hash join on the cell string, or on its 64-bit integer from h3.str_to_int — produces at most one row per pair: the measured candidate table had no duplicated hotel–station combinations.
5. Refine by great-circle distance
The merge returned 45,650 candidates. Filtering on haversine distance kept 13,931 — exactly the ground truth, with zero false pairs and zero misses:
found = set(zip(pairs.id_hotel, pairs.id_station))
print(len(found & truth), len(found - truth), len(truth - found))13931 0 0
Never skip this step. Without it, 31,719 of the 45,650 rows would be pairs more than 500 m apart. Check against something independent as well: the KD-tree truth is exact, whereas a projected spatial join is only a second opinion, as Example 3 shows.
6. Move to DuckDB when the data outgrows pandas
The same logic is one SQL statement (Example 2). Over every one of the 13,464,117 GeoNames features against 76,745 stations worldwide, it returned 101,095 pairs within 500 m, identical to a KD-tree ground truth over the whole dataset.
Code examples
Example 1 — a reusable distance join on H3 keys
import math
import h3
import pandas as pd
def h3_distance_join(left, right, distance_m, res=None, lat="lat", lon="lon",
suffixes=("_left", "_right")):
"""All pairs of points within distance_m metres, found through H3 keys and refined exactly."""
if res is None: # the finest resolution whose edge is still at least a third of the distance
res = max(r for r in range(16) if h3.average_hexagon_edge_length(r, "m") >= distance_m / 3)
k = math.ceil(distance_m / (1.5 * h3.average_hexagon_edge_length(res, "m"))) + 1
l = left.assign(_cell=[h3.latlng_to_cell(y, x, res) for y, x in zip(left[lat], left[lon])])
r = right.assign(_cell=[h3.grid_disk(h3.latlng_to_cell(y, x, res), k)
for y, x in zip(right[lat], right[lon])]).explode("_cell")
m = l.merge(r, on="_cell", suffixes=suffixes)
m["distance_m"] = haversine_m(m[lat + suffixes[0]], m[lon + suffixes[0]],
m[lat + suffixes[1]], m[lon + suffixes[1]])
kept = m[m.distance_m <= distance_m].drop(columns="_cell")
print(f"res {res}, k={k}: {len(m):,} candidates -> {len(kept):,} pairs")
return kept.reset_index(drop=True)
near = h3_distance_join(hotels, stations, 500, suffixes=("_hotel", "_station"))
wider = h3_distance_join(hotels, stations, 2000, suffixes=("_hotel", "_station"))res 9, k=3: 45,650 candidates -> 13,931 pairs
res 7, k=2: 480,233 candidates -> 113,670 pairs
The 500 m join took 0.88 s and the 2 km join 0.45 s, cell assignment included. Choosing the resolution from the distance keeps k at 2 or 3, so the key count stays manageable whatever distance is asked for. Pass the smaller table as right.
Example 2 — the same join in DuckDB, at global scale
import duckdb
def duckdb_distance_join(parquet, left_where, right_where, distance_m, res=8, k=2, threads=4):
con = duckdb.connect()
con.execute(f"set threads = {threads}")
con.execute("install h3 from community; load h3")
return con.execute(f"""
with r as (
select geonameid as rid, lat as rlat, lon as rlon,
unnest(h3_grid_disk(h3_latlng_to_cell(lat, lon, {res}), {k})) as cell
from read_parquet('{parquet}') where {right_where}
),
l as (
select geonameid as lid, lat, lon, h3_latlng_to_cell(lat, lon, {res}) as cell
from read_parquet('{parquet}') where {left_where}
)
select lid, rid,
2 * 6371008.8 * asin(sqrt(
pow(sin(radians(rlat - lat) / 2), 2)
+ cos(radians(lat)) * cos(radians(rlat)) * pow(sin(radians(rlon - lon) / 2), 2)
)) as distance_m
from l join r using (cell)
where distance_m <= {distance_m}
""").df()
near_stations = duckdb_distance_join(
"geonames_all.parquet",
left_where="feature_code not in ('RSTN', 'RSTP')",
right_where="feature_code in ('RSTN', 'RSTP')",
distance_m=500,
)
print(f"{len(near_stations):,} pairs")101,095 pairs
That took 1.67 s on four threads and 0.71 s on six, reading the Parquet file both times. A KD-tree over the same 13.4 million points found the same 101,095 pairs in 7.69 s including the load. No projected CRS could have served this query: the data spans every continent.
Example 3 — comparing with a projected spatial join
import time
import geopandas as gpd
def compare_with_sjoin(left, right, distance_m, epsg):
"""Pairs a projected-CRS sjoin finds against the true great-circle pairs."""
t = time.perf_counter()
gl = gpd.GeoDataFrame(left, geometry=gpd.points_from_xy(left.lon, left.lat), crs=4326).to_crs(epsg)
gr = gpd.GeoDataFrame(right, geometry=gpd.points_from_xy(right.lon, right.lat), crs=4326).to_crs(epsg)
joined = gpd.sjoin(gl, gr, predicate="dwithin", distance=distance_m)
seconds = time.perf_counter() - t
projected = set(zip(joined.id_left, joined.id_right))
tl, tr = cKDTree(unit_xyz(left.lat, left.lon)), cKDTree(unit_xyz(right.lat, right.lon))
chord = 2 * math.sin(distance_m / 6_371_008.8 / 2)
true = {(left.id.iat[i], right.id.iat[j])
for i, js in enumerate(tl.query_ball_tree(tr, chord)) for j in js}
print(f"EPSG:{epsg} sjoin: {len(projected):,} pairs in {seconds:.2f}s; "
f"true {len(true):,}; missed {len(true - projected)}, false {len(projected - true)}")
return true - projected, projected - true
missed, false = compare_with_sjoin(hotels, stations, 500, 3035)EPSG:3035 sjoin: 13,906 pairs in 0.14s; true 13,931; missed 43, false 18
The spatial join was six times faster than the pandas H3 join, and 61 pairs wrong. Every one was borderline: the missed pairs were 494.4–500.0 m apart on the sphere and 500.1–505.5 m apart in EPSG:3035. One, near Kraków, measured 499.9 m true and 500.6 m projected.
Explanation
Why a same-cell join fails at the edges
Two points share a cell only if no cell edge runs between them. At res 9 cell centres are about 348 m apart, so two points 500 m apart almost always sit in different cells — hence 15.9% recall. A key join treats "different key" as "not near", which is exactly the assumption a grid cannot support.
The problem is not specific to H3. Any grid key has it — geohash, S2 or a rounded coordinate — because it is the edge effect of quantising location.
Why the rule has the shape it does
A point can sit anywhere inside its cell, up to about one edge length from the centre. So a disk around its cell reaches the join distance in every direction only if its narrowest radius covers the distance plus that offset.
The narrowest radius of a disk grows steadily with k: at res 8, 969 m for k = 1, 2,421 m for k = 3, 4,598 m for k = 6 — 726 m per ring, or 1.37 edge lengths where an edge is 531 m. The rule divides the distance by 1.5 edges, rounds up, and adds a ring for the offset. It is slightly conservative — at res 7 and 1 km, k − 1 also missed nothing — which is the right direction to err.
Why a projected spatial join is not the exact answer
sjoin(predicate="dwithin") measures distance in the CRS's units. EPSG:3035 is an equal-area projection for Europe, and equal area does not mean true distance: scale varies with direction and position. The error is a few metres at 500 m, and any pair within a few metres of the threshold can flip.
On 977,803 European populated places against the same stations the pattern held: the projected sjoin took 0.82 s and missed 47 true pairs while adding 13; the H3 disk join took 0.64–0.81 s at res 8 or 9 and missed at most one pair, even with one ring fewer than the rule asks for. For most analyses 0.4% at the margin is harmless. For a regulatory buffer it is not, and for data spanning the globe there is no single CRS to choose.
Why speed is not the argument
GeoPandas' spatial index made the projected join fast — 0.14 s against 0.88 s. The case for the H3 join is elsewhere: it is exact on the sphere, it works anywhere on Earth without a CRS, the keys can be computed once and stored, and the join itself is a plain equality that every engine, from pandas to DuckDB to a data warehouse, already optimises.
Edge cases or notes
- k = 1 is often nearly right, and that is the danger. At res 8 and 500 m it missed 2 of 13,931 pairs — rare enough to pass a casual check.
- Expand the smaller side. Keys scale with rows × (1 + 3k(k + 1)).
- Pentagons need no special case. Their disks hold a few cells fewer, and 24,000 test pairs around them lost nothing with the rule's k.
- Identical coordinates are pairs at 0 m. GeoNames has many duplicated points, so decide whether a point may pair with itself before joining a table to itself.
- Nearest-neighbour questions are different. A distance join returns every pair within d; for the single nearest, grow k until a candidate appears, then search again with the k the rule gives for that candidate's distance.
- One resolution on both sides. A res-8 cell never equals a res-9 cell, and the merge will silently return nothing.
- The rule uses the average edge. It held everywhere tested, including high latitudes and pentagons, because the extra ring absorbs H3's modest variation in cell size.
Internal links
- How to perform a spatial join in Python — the geometry-based alternative
- Nearest-neighbour joins explained — distance, ties and search radius
- How to use H3 neighbours for k-ring smoothing and buffers — grid_disk, and why a disk is not a circle
- How to assign points to H3 cells in Python — computing and storing the keys
- How to use H3 in DuckDB for grid aggregation at scale — the extension used in Example 2
- How to choose the right projected CRS — when a projected spatial join is good enough
- H3 hierarchy explained — why keys from different resolutions never match
FAQ
Can I join two point datasets just by merging on the H3 cell?
Only if you want pairs that happen to share a cell. For pairs within a distance, a same-cell join found 15.9% of true hotel–station pairs within 500 m at res 9 and 59.6% at res 8.
How do I choose k for a distance?
Use k = ceil(d / (1.5 × edge)) + 1 with the resolution's average edge length. It missed no pairs in 140,000 adversarial tests, while one ring fewer missed up to 75 of 20,000.
Do I still need to compute distances after the H3 join?
Yes. The disk returns candidates in a rough hexagon; at res 9 and k = 3 only 13,931 of 45,650 candidates were within 500 m.
Is the H3 join faster than a GeoPandas spatial join?
Not necessarily. On European hotels and stations the projected sjoin took 0.14 s against 0.88 s. The H3 join's advantages are exactness on the sphere, no CRS, and portability to any engine.
Does this work for worldwide data?
Yes, which is its strongest case. In DuckDB, 13.46 million features against 76,745 stations returned 101,095 pairs in 0.71 s, identical to an exact KD-tree result.