Geospatial cell index (geohash, S2, H3)

Deep dive

How it works

In: (lat, lon) in degrees and a level L; for queries also a box or radius · Out: a 64-bit cell id; for a query, a sorted list of id ranges · State: none · Cost: point to id O(L), a few ns; a query is k range scans

1. Scale lon over 360 and lat over 180 by 2^L: integer cell coordinates (x, y) at level L.
2. Interleave the bits of x and y, lon bit first (Z-order, the geohash order). A longer id extends its ancestor's bits. S2 uses the same quadtree on a cube face but orders it along a Hilbert curve; H3 uses base-7 digits on an icosahedron. See space-filling curves and quadtrees.
3. Encode the level in the id. S2 appends a 1 bit after the path, so the level is the trailing-zero count and one integer carries it; the code copies that layout. Geohash instead relies on string length.
4. Descendants of cell c with lowest set bit lsb occupy exactly [c-(lsb-1), c+(lsb-1)].
5. Box query: choose a level where the box spans a few cells, list the cells, sort their ranges, merge adjacent ones, scan each range, then filter with the exact coordinate test.

level 2, ids by (x,y)       cells 3 and 9 touch in space, 6 ids apart
y3  5  7 13 15
y2  4  6 12 14
y1  1  3  9 11
y0  0  2  8 10
   x0 x1 x2 x3

Code

Plain Rust: 60 payload bits plus a sentinel. A leaf id is exactly a 12-char geohash (tested against u4pruydqqvj and ezs42).

// Quadtree cell ids on a lat/lon grid: Z-order bits, S2-style layout (sentinel bit).
const MAXL: u32 = 30;

fn spread(v: u32) -> u64 {
    let mut x = v as u64;
    x = (x | x << 16) & 0x0000_FFFF_0000_FFFF;
    x = (x | x << 8) & 0x00FF_00FF_00FF_00FF;
    x = (x | x << 4) & 0x0F0F_0F0F_0F0F_0F0F;
    x = (x | x << 2) & 0x3333_3333_3333_3333;
    (x | x << 1) & 0x5555_5555_5555_5555
}
fn compact(mut x: u64) -> u32 {
    x &= 0x5555_5555_5555_5555;
    x = (x | x >> 1) & 0x3333_3333_3333_3333;
    x = (x | x >> 2) & 0x0F0F_0F0F_0F0F_0F0F;
    x = (x | x >> 4) & 0x00FF_00FF_00FF_00FF;
    x = (x | x >> 8) & 0x0000_FFFF_0000_FFFF;
    ((x | x >> 16) & 0xFFFF_FFFF) as u32
}
/// Cell at `level` from grid indexes (x = lon, y = lat; `level` bits each).
fn make(level: u32, x: u32, y: u32) -> u64 {
    let z = spread(x) << 1 | spread(y); // lon bit is the higher one, like geohash
    (z << (61 - 2 * level)) | (1u64 << (60 - 2 * level)) // payload, then sentinel 1
}
fn grid(v: f64, lo: f64, span: f64, level: u32) -> u32 {
    let n = 1u64 << level;
    (((v - lo) / span * n as f64) as u64).min(n - 1) as u32
}
pub fn leaf(lat: f64, lon: f64) -> u64 {
    make(MAXL, grid(lon, -180.0, 360.0, MAXL), grid(lat, -90.0, 180.0, MAXL))
}
pub fn level(id: u64) -> u32 { MAXL - id.trailing_zeros() / 2 }
fn xy(id: u64) -> (u32, u32, u32) {
    let l = level(id);
    let z = id >> (61 - 2 * l);
    (l, compact(z >> 1), compact(z))
}
pub fn parent(id: u64, level: u32) -> u64 {
    let lsb = 1u64 << (60 - 2 * level);
    (id & lsb.wrapping_neg()) | lsb
}
/// All descendants (and the cell itself) have ids in this inclusive range.
pub fn range(id: u64) -> (u64, u64) {
    let lsb = id & id.wrapping_neg();
    (id - (lsb - 1), id + (lsb - 1))
}
/// (lat0, lat1, lon0, lon1)
pub fn bounds(id: u64) -> (f64, f64, f64, f64) {
    let (l, x, y) = xy(id);
    let (dx, dy) = (360.0 / (1u64 << l) as f64, 180.0 / (1u64 << l) as f64);
    (-90.0 + y as f64 * dy, -90.0 + (y + 1) as f64 * dy, -180.0 + x as f64 * dx, -180.0 + (x + 1) as f64 * dx)
}
/// Up to 8 same-level neighbours: longitude wraps, latitude stops at the poles.
pub fn neighbors(id: u64) -> Vec<u64> {
    let (l, x, y) = xy(id);
    let n = 1i64 << l;
    let mut out = Vec::new();
    for (dx, dy) in [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] {
        let (nx, ny) = ((x as i64 + dx).rem_euclid(n), y as i64 + dy);
        if (0..n).contains(&ny) {
            let c = make(l, nx as u32, ny as u32);
            if c != id && !out.contains(&c) { out.push(c); }
        }
    }
    out
}
/// Sorted, merged id ranges covering a lat/lon box (no dateline crossing) with level-`level` cells.
pub fn cover(lat0: f64, lat1: f64, lon0: f64, lon1: f64, level: u32) -> Vec<(u64, u64)> {
    let (x0, x1) = (grid(lon0, -180.0, 360.0, level), grid(lon1, -180.0, 360.0, level));
    let (y0, y1) = (grid(lat0, -90.0, 180.0, level), grid(lat1, -90.0, 180.0, level));
    let mut r: Vec<(u64, u64)> = (x0..=x1).flat_map(|x| (y0..=y1).map(move |y| range(make(level, x, y)))).collect();
    r.sort_unstable();
    let mut m: Vec<(u64, u64)> = Vec::new();
    for (a, b) in r {
        match m.last_mut() {
            Some(last) if a <= last.1 + 2 => last.1 = last.1.max(b), // leaf ids are odd: +2 is adjacent
            _ => m.push((a, b)),
        }
    }
    m
}
/// The leaf's 60 payload bits are exactly a 12-char geohash; c chars = top 5c bits.
pub fn geohash(leaf_id: u64, chars: u32) -> String {
    const A: &[u8; 32] = b"0123456789bcdefghjkmnpqrstuvwxyz";
    let z = leaf_id >> 1;
    (0..chars).map(|i| A[((z >> (55 - 5 * i)) & 31) as usize] as char).collect()
}

Tested: parent/range containment, 8 neighbours with dateline wrap and pole clamp, and 200k random points where box scans at levels 4/6/8 return exactly the brute-force hits.

Complexity and constants

• Encode: two bit-spreads and shifts, branch-free; order of 3 ns per point measured here (Apple M-series, release build).
• Storage: 8 bytes per id against up to 12 bytes for a geohash string; sorted neighbours share high bits and compress well.
• Query cost = ranges (seeks) + candidates (sequential reads) + exact tests. Measured, 200k uniform points, a 5 x 7 degree box with 113 hits: level 4 gave 2 ranges and 1581 candidates; level 6, 5 and 458; level 8, 14 and 163.
• Cell counts: 4^L here; 6 * 4^30 S2 leaves; H3 res 0 has 122 cells, about 7x more per resolution.

Efficiency

• Cost: point to id is O(L) bit work, a few ns without trigonometry (S2 states about 0.15 us per unit vector to id, an old measurement). A query is k range seeks plus candidate reads and exact tests. Bound: seeks at coarse levels, candidate reads and tests at fine levels.
• Levers:
  • merge adjacent ranges — one scan instead of one per cell (done in cover) — any covering
  • cap the cell count (S2 max_cells) — fewer seeks for some false positives, which are cheap sequential reads — store with costly seeks
  • Hilbert instead of Z order — fewer ranges per region; S2 states close ids stay close, helping locality — region queries
  • pre-aggregate by one fixed-level prefix (GROUP BY id >> k) — skips covering at query time — fixed-size heatmaps and joins
• Naive → standard → state of the art:
  • scan all points and run haversine on each · O(n) per query
  • Z-order cell ids, merged range scans, exact re-test · k seeks
  • S2 or H3 covering with a cell cap, pre-aggregated at a fixed level · seeks bounded by the cap

Implementation notes

• Index the 8-byte leaf id as the key; any ancestor is (id & -lsb) | lsb with the new lsb, no extra columns.
• Pick the cover level from the box size (about 4 to 16 cells across), not a constant.
• Replace spread/compact with BMI2 pdep/pext on x86 (core::arch::x86_64::_pdep_u64) if profiling shows it matters.
• Keep lat/lon next to the id so the final filter needs no join.
• This layout uses 61 bits and fits i64; S2 uses all 64, so signed columns need a cast that preserves order.

Pitfalls

• Close points do not always share a prefix: neighbours across a high-level boundary share almost nothing. Cover the box or the 8 neighbours, never one prefix.
• Geohash cells are degree rectangles: width shrinks with cos(lat). Size radius searches by the narrower side.
• Boxes crossing lon 180 need two covers; latitude has no wrap, pole cells need clamping.
• Cells are only roughly equal: S2 level 12 spans 3.04-6.38 km2, H3 hexagons vary in area too.
• H3 containment is approximate: a parent only approximately contains its seven children, while neighbours and cell borders at one resolution are exact. Pentagons have 5 neighbours; code assuming 6 breaks.
• Odd geohash lengths give non-square cells: 5 chars = 25 bits, 13 for lon and 12 for lat.

Variants and alternatives

In production

• Search and key-value stores: Elasticsearch geohash_grid aggregation, precision 1 to 12; Redis GEO keeps members in a sorted set.
• Data warehouses: BigQuery has S2_CELLIDFROMPOINT and S2_COVERINGCELLIDS (docs); Snowflake has H3_LATLNG_TO_CELL (docs); Databricks ships H3 SQL functions (docs).
• Marketplaces: H3 was developed for Uber's data science needs (docs).
• S2 is Google's library; H3 has bindings in many languages.

Implementations

Sources

• S2 Geometry, S2 Cells and S2 Cell Statistics. cells · statistics · code
• Uber, H3 documentation: index structure · cell statistics · indexing and containment · code
• Gustavo Niemeyer, Geohash, 2008. overview
• Elastic, Geohash grid aggregation. docs