Discrete Global Grids Explained: H3, S2, Geohash and Why Cells Beat Coordinates

Problem statement

A coordinate pair makes a poor key. Two GPS fixes at the same café differ in the sixth decimal place, so they never compare equal. Group by them and you get one group per row; join on them and nothing matches. So every "how many events per area" question turns into a geometry job: build polygons, pick a CRS, run a spatial join.

A discrete global grid swaps the coordinate for the identifier of the cell that contains it. That identifier is a short string or a 64-bit integer. It can be grouped, sorted, hashed, stored in any database and joined with an ordinary equality. There is no CRS to choose, because the grid covers the whole Earth.

Four grids are in common use: H3, S2, geohash and Web Mercator quadkeys. They are not interchangeable. Measured on the same data:

• Area. Quadkey tiles at zoom 12 cover 4.0× less ground at 60° than at the equator, and 131.7× less at 85°. H3 hexagons at resolution 5 vary 1.98×, and the variation does not follow latitude.
• Edges. Every grid splits close neighbours across cell boundaries. Of 33,048 pairs of real places less than 100 m apart, 52.3% landed in different geohash-7 cells and 21.4% in different H3 resolution-9 cells.

Choosing a grid means choosing which of these properties you can live with.

Quick answer

A grid turns a coordinate into a key. Each system does it in one call:

import h3, s2sphere, pygeohash, mercantile

lat, lng = 51.5007, -0.1246                      # Westminster
print("H3 res 9     ", h3.latlng_to_cell(lat, lng, 9))
print("S2 level 13  ", s2sphere.CellId.from_lat_lng(
    s2sphere.LatLng.from_degrees(lat, lng)).parent(13).to_token())
print("geohash 7    ", pygeohash.encode(lat, lng, precision=7))
print("quadkey z15  ", mercantile.quadkey(mercantile.tile(lng, lat, 15)))

H3 res 9      89194ad14c3ffff
S2 level 13   487604c4
geohash 7     gcpuvpm
quadkey z15   031313131130102

Once every row has a key, aggregation is a plain groupby. Measured on 236,014 GeoNames places at H3 resolution 5, indexing plus counting took 0.24 s. A spatial join of the same points against drawn hexagons took 0.47 s, not counting the time to build the hexagons. It also disagreed with the index on 21 points: 18 matched a neighbouring hexagon and 3 matched none. Every one of the 18 was within 1.7 m of a drawn edge.

Cell keys as rounded coordinates, and how four grids differ

1. Treat a cell key as a coordinate rounded onto the sphere

floor(lon / 0.1) is already a grid: it rounds a coordinate onto a lattice and the rounded value becomes the group key. A discrete global grid does the same rounding on the sphere, then gives the result a compact name.

Every system has levels. A coarser level means bigger cells and a shorter or smaller key, so each system also has a parent-and-child relationship. How exactly those levels nest varies by system, and that difference matters later.

2. Know how each of the four is built

H3 and S2 project a polyhedron onto the sphere, so their cells stay roughly the same size everywhere. Geohash and quadkeys subdivide a flat map, which makes the arithmetic trivial but shrinks the cells towards the poles. Quadkeys stop entirely at 85.05° north and south, where Web Mercator ends.

3. Compare how much the cell area varies

Measured within one level of each system:

system    level      area within the level
H3        res 5      hexagons 153.8–305.1 km² (1.98×); 12 pentagons at 127.8 km²
S2        level 7    all 98,304 cells, 3,175–6,529 km² (2.06×)
S2        level 10   100,000-cell sample, 48.8–102.0 km² (2.09×)
geohash   5 chars    23.9 km² at the equator, 11.9 at 60° (2.0×), 2.1 at 85°
quadkey   zoom 12    95.5 km² at the equator, 23.9 at 60° (4.0×), 0.7 at 85°

The spread alone is only half the story. The other half is where the variation sits. Geohash and quadkey cells shrink with latitude, so a map of counts per cell has a latitude gradient built into it. H3 and S2 cells vary with their position on a polyhedron face. By latitude band, the mean H3 resolution-5 area only moves between 247 and 265 km².

4. Compare the neighbours

A hexagon has one kind of neighbour. A square has two: four that share an edge and four that share only a corner.

H3 res 7, London       6 neighbours, centres 2.243–2.426 km away   ratio 1.08
geohash 6, London      8 neighbours, centres 0.611–0.975 km away   ratio 1.60

That 1.60 ratio has a practical cost. "The cell and its neighbours" covers a different distance in different directions, so a smoothing or buffer built from square neighbours is not round.

5. Expect the edge problem in every grid

However small the cells, two points close together can sit on opposite sides of a boundary. Measured on 327,277 places in Great Britain, France and Germany, which give 33,048 pairs under 100 m and 370,337 under 500 m (Example 4):

                          pairs < 100 m apart    pairs < 500 m apart
                          split across an edge   split across an edge
H3 res 9 (edge 201 m)          21.4%                  81.8%
H3 res 8 (edge 531 m)           9.4%                  42.1%
geohash 7 (153 m tall)         52.3%                  95.3%
geohash 6 (611 m tall)         10.2%                  48.7%
quadkey z15                     8.8%                  43.1%

With geohash it gets worse, because a shared prefix looks like closeness and is not. Two points 13.2 m apart on either side of the Greenwich meridian encode as gcpuzgrb and u10hb520, which share no characters at all.

6. Group and join on the key, not the geometry

import h3
import pandas as pd

cities["cell"] = [h3.latlng_to_cell(y, x, 5)
                  for y, x in zip(cities.lat.tolist(), cities.lon.tolist())]
counts = cities.groupby("cell").size()

Nothing here needs a CRS, a geometry column or a spatial index. The key is a string, so the same operation works in pandas, DuckDB, BigQuery, Spark or a key-value store. Two datasets keyed the same way can be joined with merge.

7. Pick one system per project and record the level

• H3 for analysis: aggregation, smoothing and neighbourhoods, where near-uniform hexagons matter.
• S2 for covering regions with cells of mixed sizes, and wherever the Google or BigQuery ecosystem already uses it.
• Geohash when the store can only index strings and prefix scans are the query.
• Quadkeys when the output really is Web Mercator tiles.

Keys from different systems, or from different levels of the same system, do not compare equal. Store the level beside the key, or put it in the column name.

Keys, areas and neighbours in every system

The data used on this page

Two extracts of GeoNames: cities500, every populated place of 500 people or more, as cities; and every populated place, building and site (feature classes P and S) in Great Britain, France and Germany from the allCountries dump, as places. The whole dump is also kept as points for one edge-case count. GeoNames is edited daily, so counts drift a little between downloads.

import os
import zipfile

import pandas as pd
import requests


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"):
    # GeoNames cities500 as saved on 8 October 2026 (rebuilt nightly); live: https://download.geonames.org/export/dump/cities500.zip
    fetch("https://data.spatialworkflow.io/geonames-cities500-2026-10-08.zip", "cities500.zip")
    with zipfile.ZipFile("cities500.zip") as z, z.open("cities500.txt") as f:
        pd.read_csv(f, sep="\t", header=None, usecols=[1, 4, 5, 8, 14], quoting=3,
                    names=["name", "lat", "lon", "country", "population"],
                    keep_default_na=False, na_values={"population": [""]}).to_parquet("cities500.parquet")
if not os.path.exists("places_gb_fr_de.parquet"):   # every GeoNames feature, about 400 MB zipped
    fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")
    with zipfile.ZipFile("allCountries.zip") as z, z.open("allCountries.txt") as f:
        every = pd.read_csv(f, sep="\t", header=None, usecols=[4, 5, 6, 8], quoting=3,
                            names=["lat", "lon", "fclass", "country"], keep_default_na=False)
    every[every.country.isin(["GB", "FR", "DE"]) & every.fclass.isin(["P", "S"])] \
        .reset_index(drop=True).to_parquet("places_gb_fr_de.parquet")
    every[["lat", "lon"]].to_parquet("geonames_all.parquet")

cities = pd.read_parquet("cities500.parquet")
places = pd.read_parquet("places_gb_fr_de.parquet")
print(f"cities: {len(cities):,}   places in GB, FR and DE: {len(places):,}")

cities: 236,014   places in GB, FR and DE: 327,277

Example 1 — every system's key for a DataFrame

import time
import h3, s2sphere, pygeohash, mercantile
import pandas as pd

KEYS = {
    "h3": lambda y, x: h3.latlng_to_cell(y, x, 9),
    "s2": lambda y, x: s2sphere.CellId.from_lat_lng(
        s2sphere.LatLng.from_degrees(y, x)).parent(13).to_token(),
    "geohash": lambda y, x: pygeohash.encode(y, x, precision=7),
    "quadkey": lambda y, x: mercantile.quadkey(mercantile.tile(x, y, 15)),
}


def add_grid_keys(df, lat="lat", lon="lon", systems=KEYS):
    """Add one key column per grid system, all from the same coordinates."""
    out = df.copy()
    ys, xs = out[lat].tolist(), out[lon].tolist()
    for name, key in systems.items():
        started = time.perf_counter()
        out[name] = [key(y, x) for y, x in zip(ys, xs)]
        print(f"{name:8} {out[name].nunique():>8,} cells  "
              f"{time.perf_counter() - started:5.2f} s")
    return out


keyed = add_grid_keys(cities)
print()
uk = keyed[(keyed.country == "GB") & keyed.name.isin(["London", "Cardiff", "Edinburgh"])]
print(uk[["name", "h3", "s2", "geohash", "quadkey"]].to_string(index=False))

h3        235,429 cells   0.34 s
s2        231,458 cells   1.41 s
geohash   235,779 cells   0.28 s
quadkey   232,202 cells   0.84 s

     name              h3       s2 geohash         quadkey
   London 89195da49b7ffff 487604cc gcpvj0u 031313131130100
Edinburgh 8919727650fffff 4887c794 gcvwr3b 031133233211123
  Cardiff 8919586b4bbffff 486e1cb4 gcjszes 031313031031310

The cell counts are for every place in GeoNames cities500.

The speed gap comes from the implementations: H3's core is C, while s2sphere is pure Python. It tells you nothing about the grids themselves.

Example 2 — the ground area of a cell in each system

import math
import h3, s2sphere, pygeohash, mercantile

R = 6371.007180918475          # km — the sphere h3.cell_area uses


def rect_km2(south, north, west, east):
    """Exact spherical area of a latitude–longitude rectangle."""
    return (R ** 2 * math.radians(east - west)
            * abs(math.sin(math.radians(north)) - math.sin(math.radians(south))))


def cell_areas_km2(lat, lon, h3_res=5, s2_level=10, gh_precision=5, zoom=12):
    """Ground area of the cell containing one point, in each system."""
    cell = h3.latlng_to_cell(lat, lon, h3_res)
    s2_cell = s2sphere.Cell(s2sphere.CellId.from_lat_lng(
        s2sphere.LatLng.from_degrees(lat, lon)).parent(s2_level))
    gy, gx, dy, dx = pygeohash.decode_exactly(
        pygeohash.encode(lat, lon, precision=gh_precision))
    tile = mercantile.bounds(mercantile.tile(lon, lat, zoom))
    return {
        "h3": h3.cell_area(cell, "km^2"),
        "s2": s2_cell.exact_area() * R ** 2,
        "geohash": rect_km2(gy - dy, gy + dy, gx - dx, gx + dx),
        "quadkey": rect_km2(tile.south, tile.north, tile.west, tile.east),
    }


print(" lat    h3 r5   s2 l10     gh 5      z12")
for lat in (0.5, 30, 45, 60, 75):
    a = cell_areas_km2(lat, 10.3)
    print(f"{lat:4}  {a['h3']:7.1f}  {a['s2']:7.1f}  {a['geohash']:7.1f}  {a['quadkey']:7.1f}")

 lat    h3 r5   s2 l10     gh 5      z12
 0.5    237.5     82.5     23.9     95.5
  30    287.6     87.9     20.7     71.6
  45    255.0     60.0     16.9     47.8
  60    180.1     84.1     11.9     23.9
  75    211.3     88.8      6.2      6.4

Along the meridian at 10.3° E, the last two columns fall steadily with latitude. The first two go up and down: the 60° hexagon is smaller than the one at 75°. For a density map this matters. Divide by the area of each cell, never by a single nominal area per level.

Example 3 — how often near neighbours are split

import math
import numpy as np
from scipy.spatial import cKDTree

EARTH_M = 6371008.8


def split_rate(lat, lon, key, distance_m):
    """Share of point pairs closer than distance_m that land in different cells."""
    lat, lon = np.asarray(lat), np.asarray(lon)
    phi, lam = np.radians(lat), np.radians(lon)
    unit = np.column_stack([np.cos(phi) * np.cos(lam),
                            np.cos(phi) * np.sin(lam), np.sin(phi)])
    chord = 2 * math.sin(distance_m / EARTH_M / 2)
    pairs = cKDTree(unit).query_pairs(chord, output_type="ndarray")
    a, b = pairs[:, 0], pairs[:, 1]
    moved = (lat[a] != lat[b]) | (lon[a] != lon[b])      # ignore duplicates
    a, b = a[moved], b[moved]
    used = np.unique(np.concatenate([a, b]))
    keys = {i: key(lat[i], lon[i]) for i in used.tolist()}
    split = sum(keys[i] != keys[j] for i, j in zip(a.tolist(), b.tolist()))
    return len(a), split / len(a)


for label, key in [("H3 res 9", KEYS["h3"]), ("geohash 7", KEYS["geohash"])]:
    n, rate = split_rate(places.lat, places.lon, key, 100)
    print(f"{label:10} {n:,} pairs within 100 m, {rate:.1%} split across a cell edge")

H3 res 9   33,048 pairs within 100 m, 21.4% split across a cell edge
geohash 7  33,048 pairs within 100 m, 52.3% split across a cell edge

Run it on your own data at the distance your analysis cares about, before you trust a same-cell join. If the split rate is not close to zero, the join also has to look in neighbouring cells.

Example 4 — the measurements quoted on this page

Every figure in the steps above is printed here, from the data loaded at the top of this section.

import sys
import numpy as np
import geopandas as gpd
from shapely.geometry import Point, Polygon

# Indexing and counting against a spatial join on drawn hexagons, H3 resolution 5
ys, xs = cities.lat.tolist(), cities.lon.tolist()
started = time.perf_counter()
cities["cell"] = [h3.latlng_to_cell(y, x, 5) for y, x in zip(ys, xs)]
counts = cities.groupby("cell").size()
t_index = time.perf_counter() - started
hexes = gpd.GeoDataFrame(
    {"drawn": counts.index},
    geometry=[Polygon([(x, y) for y, x in h3.cell_to_boundary(c)]) for c in counts.index],
    crs=4326)
pts = gpd.GeoDataFrame(cities[["cell"]], geometry=gpd.points_from_xy(cities.lon, cities.lat), crs=4326)
started = time.perf_counter()
joined = gpd.sjoin(pts, hexes, how="left", predicate="within")
t_join = time.perf_counter() - started
joined = joined[~joined.index.duplicated()]
other = joined[joined.drawn.notna() & (joined.drawn != joined.cell)]
none = joined[joined.drawn.isna()]


def metres_to_outline(point, cell):
    """Distance from a point to a cell's drawn outline, in a local equirectangular frame."""
    k = math.cos(math.radians(point.y))
    ring = Polygon([(x * k, y) for y, x in h3.cell_to_boundary(cell)]).exterior
    return ring.distance(Point(point.x * k, point.y)) * 111_195


print(f"{len(cities):,} places: index + count {t_index:.2f} s, spatial join {t_join:.2f} s")
print(f"disagree: {len(other) + len(none)} ({len(other)} in a neighbour's outline, {len(none)} in none); "
      f"neighbour cases at most {max(metres_to_outline(p, c) for p, c in zip(other.geometry, other.drawn)):.1f} m from the edge")

# How far a straight lon/lat edge sits from H3's true (great-circle) edge, resolution 5
def unit(lat, lng):
    a, b = math.radians(lat), math.radians(lng)
    return np.array([math.cos(a) * math.cos(b), math.cos(a) * math.sin(b), math.sin(a)])


gaps = []
for c in counts.index[::10]:
    ring = h3.cell_to_boundary(c)
    for (y1, x1), (y2, x2) in zip(ring, ring[1:] + ring[:1]):
        if abs(x1 - x2) > 180:
            continue
        true_mid = unit(y1, x1) + unit(y2, x2)
        true_mid /= np.linalg.norm(true_mid)
        flat_mid = unit((y1 + y2) / 2, (x1 + x2) / 2)
        gaps.append(2 * math.asin(min(1, np.linalg.norm(true_mid - flat_mid) / 2)) * EARTH_M)
print(f"edge midpoints, drawn against true: median {np.median(gaps):.1f} m, max {max(gaps):.1f} m "
      f"over {len(gaps):,} edges")

# Step 3: area within one level of each system
res5 = [c for base in h3.get_res0_cells() for c in h3.cell_to_children(base, 5)]
area = np.array([h3.cell_area(c, "km^2") for c in res5])
pent = np.array([h3.is_pentagon(c) for c in res5])
print(f"\nH3 res 5: {len(res5):,} cells, hexagons {area[~pent].min():.1f}–{area[~pent].max():.1f} km² "
      f"({area[~pent].max() / area[~pent].min():.2f}×), {pent.sum()} pentagons at {area[pent].mean():.1f} km²")
band = np.floor(np.abs([h3.cell_to_latlng(c)[0] for c in res5]) / 15)
means = pd.Series(area).groupby(band).mean()
print(f"H3 res 5 mean area by 15° band of latitude: {means.min():.0f}–{means.max():.0f} km²")

s2_areas, cell_id = [], s2sphere.CellId.begin(7)
while cell_id != s2sphere.CellId.end(7):
    s2_areas.append(s2sphere.Cell(cell_id).exact_area() * R ** 2)
    cell_id = cell_id.next()
print(f"S2 level 7: {len(s2_areas):,} cells, {min(s2_areas):,.0f}–{max(s2_areas):,.0f} km² "
      f"({max(s2_areas) / min(s2_areas):.2f}×)")
rng = np.random.default_rng(0)
v = rng.normal(size=(100_000, 3))
v /= np.linalg.norm(v, axis=1)[:, None]
lat_r, lon_r = np.degrees(np.arcsin(v[:, 2])), np.degrees(np.arctan2(v[:, 1], v[:, 0]))
s2_10 = [s2sphere.Cell(s2sphere.CellId.from_lat_lng(s2sphere.LatLng.from_degrees(y, x)).parent(10))
         .exact_area() * R ** 2 for y, x in zip(lat_r, lon_r)]
print(f"S2 level 10: 100,000 random cells, {min(s2_10):.1f}–{max(s2_10):.1f} km² "
      f"({max(s2_10) / min(s2_10):.2f}×)")
at = {lat: cell_areas_km2(lat, 10.3) for lat in (0.01, 60, 85)}
for name, label in [("geohash", "geohash 5"), ("quadkey", "quadkey z12")]:
    eq, a60, a85 = (at[lat][name] for lat in (0.01, 60, 85))
    print(f"{label}: {eq:.1f} km² at the equator, {a60:.1f} at 60° ({eq / a60:.1f}×), "
          f"{a85:.1f} at 85° ({eq / a85:.1f}×)")

# Step 4: neighbours of a cell in London
lat, lng = 51.5074, -0.1278
home = h3.latlng_to_cell(lat, lng, 7)
d = [h3.great_circle_distance(h3.cell_to_latlng(home), h3.cell_to_latlng(n), "km")
     for n in h3.grid_disk(home, 1) if n != home]
print(f"\nH3 res 7, London   {len(d)} neighbours, centres {min(d):.3f}–{max(d):.3f} km away   "
      f"ratio {max(d) / min(d):.2f}")
gy, gx, dy, dx = pygeohash.decode_exactly(pygeohash.encode(lat, lng, precision=6))
d = [h3.great_circle_distance((gy, gx), (gy + 2 * dy * i, gx + 2 * dx * j), "km")
     for i in (-1, 0, 1) for j in (-1, 0, 1) if (i, j) != (0, 0)]
print(f"geohash 6, London  {len(d)} neighbours, centres {min(d):.3f}–{max(d):.3f} km away   "
      f"ratio {max(d) / min(d):.2f}")

# Step 5: split rates on places in Great Britain, France and Germany
split_keys = [
    (f"H3 res 9 (edge {h3.average_hexagon_edge_length(9, 'm'):.0f} m)", KEYS["h3"]),
    (f"H3 res 8 (edge {h3.average_hexagon_edge_length(8, 'm'):.0f} m)", lambda y, x: h3.latlng_to_cell(y, x, 8)),
    (f"geohash 7 ({180 / 2 ** 17 * 111_195:.0f} m tall)", KEYS["geohash"]),
    (f"geohash 6 ({180 / 2 ** 15 * 111_195:.0f} m tall)", lambda y, x: pygeohash.encode(y, x, precision=6)),
    ("quadkey z15", KEYS["quadkey"]),
]
print(f"\n{len(places):,} places")
for label, key in split_keys:
    n1, r1 = split_rate(places.lat, places.lon, key, 100)
    n5, r5 = split_rate(places.lat, places.lon, key, 500)
    print(f"{label:24} {r1:6.1%} of {n1:,} pairs < 100 m   {r5:6.1%} of {n5:,} pairs < 500 m")
a = cell_areas_km2(50, 8, h3_res=8, gh_precision=6, zoom=15)
print(f"at 50° N: quadkey z15 {a['quadkey']:.2f} km², H3 res 8 {a['h3']:.2f} km², geohash 6 {a['geohash']:.2f} km²")

# Geohash prefixes across the meridian and the equator
for (y1, x1), (y2, x2) in [((51.4779, -0.00009), (51.4779, 0.0001)), ((-0.0001, 10.3), (0.0001, 10.3))]:
    g1, g2 = pygeohash.encode(y1, x1, precision=8), pygeohash.encode(y2, x2, precision=8)
    print(f"{g1} and {g2}: {h3.great_circle_distance((y1, x1), (y2, x2), 'm'):.1f} m apart")
gy, gx, dy, dx = pygeohash.decode_exactly(pygeohash.encode(0.01, 10.3, precision=6))
print(f"geohash 6 at the equator: {2 * dy * 111.195:.3f} km tall, "
      f"{2 * dx * 111.195 * math.cos(math.radians(gy)):.3f} km wide")

# Edge-case figures: pentagons, key storage, the antimeridian
allpts = pd.read_parquet("geonames_all.parquet", columns=["lat", "lon"])
in_pentagon = 0
for p in h3.get_pentagons(5):
    py, px = h3.cell_to_latlng(p)
    near = allpts[(abs(allpts.lat - py) < 1) & (((allpts.lon - px + 180) % 360 - 180).abs() < 3)]
    in_pentagon += sum(h3.latlng_to_cell(y, x, 5) == p for y, x in zip(near.lat, near.lon))
print(f"\n{in_pentagon} of {len(allpts):,} GeoNames points fall in a resolution-5 pentagon")
two_m = pd.Series(res5[:2_000_000])
as_text = two_m.astype("string[pyarrow]").memory_usage(deep=True, index=False)
as_int = two_m.map(h3.str_to_int).astype("uint64").memory_usage(deep=True, index=False)
print(f"2,000,000 H3 keys: {as_text / 1e6:.1f} MB as Arrow strings, {as_int / 1e6:.1f} MB as uint64")
wide = sum(np.ptp([x for _, x in h3.cell_to_boundary(c)]) > 180 for c in res5)
print(f"{wide:,} resolution-5 cells span more than 180° of longitude when drawn naively")

236,014 places: index + count 0.24 s, spatial join 0.47 s
disagree: 21 (18 in a neighbour's outline, 3 in none); neighbour cases at most 1.7 m from the edge
edge midpoints, drawn against true: median 0.9 m, max 7.8 m over 50,129 edges

H3 res 5: 2,016,842 cells, hexagons 153.8–305.1 km² (1.98×), 12 pentagons at 127.8 km²
H3 res 5 mean area by 15° band of latitude: 247–265 km²
S2 level 7: 98,304 cells, 3,175–6,529 km² (2.06×)
S2 level 10: 100,000 random cells, 48.8–102.0 km² (2.09×)
geohash 5: 23.9 km² at the equator, 11.9 at 60° (2.0×), 2.1 at 85° (11.5×)
quadkey z12: 95.5 km² at the equator, 23.9 at 60° (4.0×), 0.7 at 85° (131.7×)

H3 res 7, London   6 neighbours, centres 2.243–2.426 km away   ratio 1.08
geohash 6, London  8 neighbours, centres 0.611–0.975 km away   ratio 1.60

327,277 places
H3 res 9 (edge 201 m)     21.4% of 33,048 pairs < 100 m    81.8% of 370,337 pairs < 500 m
H3 res 8 (edge 531 m)      9.4% of 33,048 pairs < 100 m    42.1% of 370,337 pairs < 500 m
geohash 7 (153 m tall)    52.3% of 33,048 pairs < 100 m    95.3% of 370,337 pairs < 500 m
geohash 6 (611 m tall)    10.2% of 33,048 pairs < 100 m    48.7% of 370,337 pairs < 500 m
quadkey z15                8.8% of 33,048 pairs < 100 m    43.1% of 370,337 pairs < 500 m
at 50° N: quadkey z15 0.62 km², H3 res 8 0.67 km², geohash 6 0.48 km²
gcpuzgrb and u10hb520: 13.2 m apart
kpzrfrup and s0p242h0: 22.2 m apart
geohash 6 at the equator: 0.611 km tall, 1.222 km wide

31 of 13,472,317 GeoNames points fall in a resolution-5 pentagon
2,000,000 H3 keys: 46.0 MB as Arrow strings, 16.0 MB as uint64
1,547 resolution-5 cells span more than 180° of longitude when drawn naively

Keys over geometries, and the edge problem no grid avoids

Why a key is cheaper than a geometry

A spatial join builds geometries, builds an index over one side, finds candidates by bounding box and then tests each candidate exactly. An equality join on a key hashes one column and looks it up. Only the second one runs in any system that can join two tables. That is why cell keys turn up in data warehouses, stream processors and key-value stores that have no geometry type.

The key also removes the CRS question. A GeoPandas dwithin join needs a projected CRS that suits the whole dataset. A global dataset has no such CRS, but it has one H3 grid.

Why no grid has equal areas and identical cells

A sphere can be tiled with identical regular polygons only as the five Platonic solids, and the finest of those has 20 faces. Any finer grid has to give something up: equal area, equal shape, or simple arithmetic.

Geohash and quadkeys keep the arithmetic. Their keys are just bit-interleaved longitude and latitude, or tile numbers, and the cost is area that collapses towards the poles. H3 and S2 keep the shape nearly uniform and accept a bounded area spread of about 2×. Equal-area designs such as ISEA and rHEALPix exist but are far less common in Python tooling.

Why hexagons make neighbourhoods simpler

Every neighbour of a hexagon shares a full edge, and its centre is about the same distance away: a ratio of 1.08 measured in London against 1.60 for geohash. So "within k steps" traces an almost circular shape, and a smoothing kernel treats all six directions alike.

The price is the hierarchy. Seven hexagons cannot tile a larger hexagon exactly, so H3 children stick out past their parent's edge. The S2 quadtree and geohash prefixes nest perfectly.

Why the edge problem cannot be designed away

A partition has boundaries, and a pair of points closer together than the cell size will sometimes straddle one. The split rate grows with the ratio of pair distance to cell size. At H3 resolution 9 it went from 21.4% at 100 m to 81.8% at 500 m.

Cell shape matters much less than cell size. Near 50° N, three cells of about the same area split 100 m pairs at about the same rate: a quadkey z15 tile (0.62 km²) split 8.8%, an H3 resolution-8 hexagon (0.67 km²) 9.4%, and a six-character geohash (0.48 km²) 10.2%. No shape brings the rate to zero. The only real fix is to search the neighbouring cells as well, then check the exact distance.

Why a join against drawn hexagons disagrees with the index

cell_to_boundary gives the cell's vertices. A GIS then joins those vertices with straight lines in longitude and latitude, but H3's true edges are computed on the sphere. At resolution 5, the midpoints of the two versions of an edge were 0.9 m apart at the median and 7.8 m at most, over 50,129 edges.

The 21 disagreeing points out of 236,014 were exactly the ones in that sliver. 18 fell inside a neighbour's drawn outline. The other 3 fell inside the drawn outline of a cell that no point had been indexed to, so they matched nothing. The index is the authority. The drawn polygon is an approximation for display.

Quadkey limits, geohash prefixes and H3 pentagons

• Quadkeys end at 85.05°. Web Mercator does not reach the poles, so anything beyond that latitude has no tile.
• A shared geohash prefix is not proximity. Points 13.2 m apart across the Greenwich meridian share zero characters, and so do points 22.2 m apart across the equator.
• Geohash cells are not square. At six characters a cell is 0.611 km tall and 1.222 km wide at the equator. Odd and even precisions alternate between near-square and 2:1.
• H3 has 12 pentagons at every resolution. They have five neighbours and all sit in the ocean. 31 of 13.47 million GeoNames points fell inside a resolution-5 pentagon.
• Store the level with the key. An H3 index carries its resolution, but a geohash only carries its length, and a quadkey's length is its zoom.
• Integers are cheaper than strings. Two million H3 keys took 46.0 MB as Arrow strings and 16.0 MB as uint64.
• Drawn cells crossing the antimeridian wrap around the map. 1,547 H3 resolution-5 cells span more than 180° of longitude when their vertices are joined naively.
• Keys do not convert between systems. Going from a geohash to an H3 cell means going back through a coordinate, and you lose whatever precision the first key had already thrown away.

Internal links

• H3 explained: how a hexagonal grid indexes the whole Earth — the grid most Python analysis uses
• Geohash explained: prefixes, precision and the edge problem — the string grid in detail
• S2 cells explained: squares on a cube — when quadrilaterals win
• How to index data by quadkey and web map tile — the tile-aligned grid
• Choosing an H3 resolution — cell size against what the data supports
• How to join two point datasets on an H3 index — the fix for the edge problem
• How to aggregate millions of points into a grid with DuckDB — the degree-grid version of the same idea
• Spatial indexes explained: R-trees — the geometry-based alternative
• The modifiable areal unit problem explained — why the cell size changes the answer

FAQ

What is a discrete global grid system?

A way of dividing the whole Earth into cells with unique identifiers, at several levels of detail. It turns a coordinate into a key that can be grouped, joined and stored without geometry or a CRS.

Is H3 better than geohash?

For analysis, usually yes. Its hexagons vary 1.98× in area with no latitude trend, and their neighbours are almost equidistant. Geohash is simpler and supports prefix scans, but its cells lose area towards the poles and a shared prefix does not mean two points are close.

Are H3 cells equal area?

No. At resolution 5 the hexagons range from 153.8 to 305.1 km² and the pentagons are 127.8 km². Divide counts by each cell's own area when densities matter.

Why do two nearby points get different cell keys?

They are on opposite sides of a cell edge. Of pairs of places less than 100 m apart, 21.4% were split at H3 resolution 9 and 52.3% at geohash precision 7.

Can I convert a geohash to an H3 cell?

Not directly. Decode the geohash to its centre and index that point, keeping in mind that the geohash already discarded the position inside its cell.

Do I still need a CRS?

Not to build keys, group or join. You do need one when you draw the cells on a map, or measure distances and areas outside what the grid library provides.