How to Join Two Point Datasets on an H3 Index Instead of a Spatial Join
On this page
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,601 GeoNames hotels and 33,572 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.
A ground truth, then the k that makes the join exact
1. Load both sides and build a ground truth once
import math
import duckdb
import h3
import numpy as np
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()
print(f"{len(hotels):,} hotels, {len(stations):,} stations")102,601 hotels, 33,572 stations
That gives 102,601 hotels and 33,572 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)])
distance = 500
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
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,641 of 50,164 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 603 pairs (res 7 at 2 km; Example 4). 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,572 stations became 1,242,164 keys. Expanding the hotels instead would have made 3,796,237.
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,653 candidates. Filtering on haversine distance kept 13,931 — exactly the ground truth, with zero false pairs and zero misses:
import math
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 # the quick answer, on the loaded data
k = math.ceil(distance / (1.5 * h3.average_hexagon_edge_length(res, "m"))) + 1
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")
candidates = a.merge(b, on="cell", suffixes=("_hotel", "_station"))
pairs = candidates[haversine_m(candidates.lat_hotel, candidates.lon_hotel,
candidates.lat_station, candidates.lon_station) <= distance]
print(k, len(b), len(candidates), len(pairs))
found = set(zip(pairs.id_hotel, pairs.id_station))
print(len(found & truth), len(found - truth), len(truth - found))3 1242164 45653 13931
13931 0 0
Never skip this step. Without it, 31,722 of the 45,653 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,472,317 GeoNames features against 76,833 stations worldwide, it returned 101,227 pairs within 500 m, identical to a KD-tree ground truth over the whole dataset.
Hotels to stations on H3 keys, in Python and DuckDB
The data used on this page
Every feature in the GeoNames allCountries dump with its feature code: hotels are HTL, railway stations RSTN and RSTP. GeoNames is edited daily, so counts drift a little.
import os
import zipfile
import pandas as pd
import requests
if not os.path.exists("geonames_all.parquet"): # about 400 MB zipped
if not os.path.exists("allCountries.zip"):
with requests.get("https://download.geonames.org/export/dump/allCountries.zip", stream=True, timeout=600) as r:
r.raise_for_status()
with open("allCountries.zip", "wb") as out:
for chunk in r.iter_content(1 << 20):
out.write(chunk)
with zipfile.ZipFile("allCountries.zip") as z, z.open("allCountries.txt") as f:
pd.read_csv(f, sep="\t", header=None, usecols=[0, 4, 5, 7], quoting=3,
names=["geonameid", "lat", "lon", "feature_code"],
dtype={"feature_code": str}).to_parquet("geonames_all.parquet")
print(len(pd.read_parquet("geonames_all.parquet", columns=["geonameid"])), "GeoNames features")13472317 GeoNames features
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)
import time
start = time.perf_counter()
near = h3_distance_join(hotels, stations, 500, suffixes=("_hotel", "_station"))
print(f"500 m join: {time.perf_counter() - start:.2f} s")
start = time.perf_counter()
wider = h3_distance_join(hotels, stations, 2000, suffixes=("_hotel", "_station"))
print(f"2 km join: {time.perf_counter() - start:.2f} s")res 9, k=3: 45,653 candidates -> 13,931 pairs
500 m join: 0.63 s
res 7, k=2: 480,260 candidates -> 113,674 pairs
2 km join: 0.39 s
The 500 m join took 0.63 s and the 2 km join 0.39 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")
every = con.execute("select lat, lon, feature_code from read_parquet('geonames_all.parquet') "
"where lat is not null").df()
is_station = every.feature_code.isin(["RSTN", "RSTP"]).to_numpy()
is_other = (every.feature_code.notna() & ~every.feature_code.isin(["RSTN", "RSTP"])).to_numpy() # as SQL's NOT IN
tree_all = cKDTree(unit_xyz(every.lat, every.lon))
tree_st = cKDTree(unit_xyz(every.lat[is_station], every.lon[is_station]))
exact = sum(int(is_other[js].sum()) for js in tree_st.query_ball_tree(tree_all, chord))
print(f"{is_station.sum():,} stations worldwide; KD-tree ground truth {exact:,} pairs; "
f"{(~is_station & ~is_other).sum():,} features with no feature code sit on neither side")101,227 pairs
76,833 stations worldwide; KD-tree ground truth 101,227 pairs; 94,926 features with no feature code sit on neither side
A KD-tree over the same 13.5 million points found the same 101,227 pairs, once features with no feature code were left out of both — SQL's NOT IN drops them silently, and a first comparison that kept them found 128 more. 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,907 pairs in 0.16s; true 13,931; missed 43, false 19
The spatial join was several times faster than the pandas H3 join, and 62 pairs wrong. Every one was borderline: the missed pairs were 494.4–500.0 m apart on the sphere, just inside the threshold, where a projected metre is not quite a true one.
Example 4 — testing the rule for k where it is most likely to fail
from pyproj import Geod
geod = Geod(ellps="WGS84")
rng = np.random.default_rng(0)
def rule(res, distance):
return math.ceil(distance / (1.5 * h3.average_hexagon_edge_length(res, "m"))) + 1
def misses(res, distance, k, origins):
"""Pairs placed at exactly `distance` whose far point falls outside the k-disk of the near one."""
lats, lons = origins
bearing = rng.uniform(0, 360, len(lats))
lon2, lat2, _ = geod.fwd(lons, lats, bearing, np.full(len(lats), distance))
return sum(h3.latlng_to_cell(b, a, res) not in set(h3.grid_disk(h3.latlng_to_cell(y, x, res), k))
for y, x, b, a in zip(lats, lons, lat2, lon2))
total = 0
for res, distance in [(7, 2000), (8, 500), (8, 2000), (9, 500), (9, 1000), (10, 250), (10, 500)]:
lats = np.degrees(np.arcsin(rng.uniform(-0.95, 0.95, 20_000)))
lons = rng.uniform(-180, 180, 20_000)
k = rule(res, distance)
total += 20_000
print(f"res {res:2} at {distance:>5} m: k={k} missed {misses(res, distance, k, (lats, lons))}, "
f"k-1 missed {misses(res, distance, k - 1, (lats, lons))}")
pent_lats, pent_lons = zip(*[h3.cell_to_latlng(p) for p in h3.get_pentagons(8) for _ in range(2000)])
print(f"{total:,} random pairs; around the 12 res-8 pentagons, {len(pent_lats):,} pairs at 499 m: "
f"missed {misses(8, 499, rule(8, 499), (np.array(pent_lats), np.array(pent_lons)))}")res 7 at 2000 m: k=2 missed 0, k-1 missed 603
res 8 at 500 m: k=2 missed 0, k-1 missed 1
res 8 at 2000 m: k=4 missed 0, k-1 missed 58
res 9 at 500 m: k=3 missed 0, k-1 missed 42
res 9 at 1000 m: k=5 missed 0, k-1 missed 53
res 10 at 250 m: k=4 missed 0, k-1 missed 1
res 10 at 500 m: k=6 missed 0, k-1 missed 209
140,000 random pairs; around the 12 res-8 pentagons, 24,000 pairs at 499 m: missed 0
Why same-cell joins fail at the edges
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, as smoothing over H3 neighbours measured it: 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.
For most analyses a handful of pairs 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.12 s against 0.63 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.
Nearly-right k and duplicate coordinates
- 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,653 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.12 s against 0.63 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.47 million features against 76,833 stations returned 101,227 pairs, identical to an exact KD-tree result.