Geohash Bucketing for Point-in-Zone Joins

Within the Local Spatial Processing Patterns guide, this page describes an alternative to the custom grid index built in spatial joins in constrained environments: instead of hashing raw coordinates into a bespoke (cx, cy) tuple, encode each fix into a geohash cell and use that cell — plus its neighbours — as the bucket key against a dictionary of zone assignments. The appeal is not that geohash is faster than a custom grid; it is that the encoding is a well-known, self-describing scheme with a fixed precision-to-cell-size relationship, which makes it easy to reason about, easy to debug in the field with off-the-shelf tools, and trivial to serialize as a single integer instead of a coordinate pair.

The join this page targets is the common “is this asset currently inside zone X” problem: matching a stream of GPS fixes against a fixed dictionary of named zones — delivery areas, maintenance sectors, exclusion boundaries — where the zone geometry itself is small enough, or coarse enough, that a bucket-level match is often sufficient without an expensive containment predicate behind it.

Bucketing scheme selection rationale #

A geohash encodes latitude and longitude by recursively bisecting each axis and interleaving the resulting bits into a single value. Every additional bit of precision halves the cell’s extent along one axis, alternating between longitude and latitude, so cell size shrinks predictably and by design as precision increases:

The right precision is the one whose cell width is close to, but larger than, the zone geometry’s smallest meaningful feature — too coarse and a cell straddles multiple zones, forcing every lookup to disambiguate; too fine and a single zone spans thousands of cells, bloating the dictionary and multiplying the neighbour cells a boundary-adjacent fix has to check.

Three properties make this scheme fit the constraint envelope of a gateway node:

1. Integer keys, not strings. The traditional base32 geohash string ("9q8yy") is convenient for humans and URLs but costly for a hot-path dictionary key: string hashing and comparison cost more than integer hashing, and each string is a separate heap object. Truncating the interleaved bit pattern to a fixed-width integer keeps the same cell semantics while making the bucket key a single machine word.
2. Neighbour cells, not just the exact cell. A fix landing within a few meters of a cell boundary can be geometrically inside a zone whose geohash cell differs from the fix’s own cell, purely because the grid lines do not know about zone boundaries. The standard fix is to check the fix’s own cell plus its eight neighbours, trading a small constant-factor increase in lookups for eliminating boundary misses.
3. Precision chosen against GPS drift, not zone size alone. A consumer GNSS receiver drifts several meters under normal conditions and tens of meters under multipath; picking a precision whose cell width is smaller than that drift means the same physical position can hash to different cells between consecutive fixes, which is exactly why the neighbour-cell check exists rather than being an optional refinement.

Checking a fix’s cell alone misses zones that begin just across a grid line; the eight neighbours close that gap.

The complete implementation #

The encoder below interleaves latitude and longitude bits directly into a Python integer rather than producing a base32 string, and exposes the neighbour cells needed to cover boundary-adjacent fixes:

# geohash_bucket.py — integer geohash bucketing for point-in-zone joins.
# Pure stdlib; no invented APIs. Precision is expressed in bits, split evenly
# (rounding favors longitude) between the two axes, alternating per bit as in
# the standard geohash bit-interleaving scheme.

PRECISION_BITS = 30  # ~1.2 km x 0.6 km cells; tune against zone size and GPS drift

def _clamp(value: float, lo: float, hi: float) -> float:
    return lo if value < lo else hi if value > hi else value

def encode_cell(lat: float, lon: float, bits: int = PRECISION_BITS) -> int:
    """Interleave lon/lat bits into a single integer bucket key. Longitude
    takes the even bit positions, latitude the odd ones, matching the
    standard geohash interleave order."""
    lat = _clamp(lat, -90.0, 90.0)
    lon = _clamp(lon, -180.0, 180.0)
    lat_lo, lat_hi = -90.0, 90.0
    lon_lo, lon_hi = -180.0, 180.0
    cell = 0
    for i in range(bits):
        cell <<= 1
        if i % 2 == 0:
            mid = (lon_lo + lon_hi) / 2.0
            if lon >= mid:
                cell |= 1
                lon_lo = mid
            else:
                lon_hi = mid
        else:
            mid = (lat_lo + lat_hi) / 2.0
            if lat >= mid:
                cell |= 1
                lat_lo = mid
            else:
                lat_hi = mid
    return cell

def cell_bounds(cell: int, bits: int = PRECISION_BITS) -> tuple[float, float, float, float]:
    """Reverse encode_cell to recover the cell's (min_lon, min_lat, max_lon, max_lat)."""
    lat_lo, lat_hi = -90.0, 90.0
    lon_lo, lon_hi = -180.0, 180.0
    for i in range(bits - 1, -1, -1):
        bit = (cell >> i) & 1
        if (bits - 1 - i) % 2 == 0:
            mid = (lon_lo + lon_hi) / 2.0
            lon_lo, lon_hi = (mid, lon_hi) if bit else (lon_lo, mid)
        else:
            mid = (lat_lo + lat_hi) / 2.0
            lat_lo, lat_hi = (mid, lat_hi) if bit else (lat_lo, mid)
    return lon_lo, lat_lo, lon_hi, lat_hi

def neighbour_cells(lat: float, lon: float, bits: int = PRECISION_BITS) -> set[int]:
    """Return the fix's own cell plus the eight surrounding cells, computed by
    re-encoding points offset by one cell width in each direction rather than
    walking bit patterns -- simpler to get right and to audit."""
    min_lon, min_lat, max_lon, max_lat = cell_bounds(encode_cell(lat, lon, bits), bits)
    width, height = max_lon - min_lon, max_lat - min_lat
    keys = set()
    for dlon in (-width, 0.0, width):
        for dlat in (-height, 0.0, height):
            probe_lat = _clamp(lat + dlat, -90.0, 90.0)
            probe_lon = ((lon + dlon + 180.0) % 360.0) - 180.0  # wrap the antimeridian
            keys.add(encode_cell(probe_lat, probe_lon, bits))
    return keys

def lookup_zones(zone_dict: dict[int, list[str]], lat: float, lon: float,
                  bits: int = PRECISION_BITS) -> set[str]:
    """Hash-lookup join: union the zone ids registered against the fix's cell
    and its neighbours. A cheap bucket-level result -- callers needing exact
    containment should confirm hits with the predicate from spatial joins in
    constrained environments before acting on a match."""
    matched: set[str] = set()
    for key in neighbour_cells(lat, lon, bits):
        matched.update(zone_dict.get(key, ()))
    return matched

neighbour_cells recomputes cell membership by offsetting the probe point, not by manipulating the bit pattern directly. That trade costs a handful of extra float comparisons per fix but makes the antimeridian wrap and pole clamping trivially correct, which is worth far more on a field device than shaving microseconds off a lookup that already runs in well under a millisecond.

Constraint validation table #

Gotchas and edge cases #

• GPS drift can move a fix across a cell boundary between consecutive readings even while the device is stationary. This is not a bug in the encoder; it is why the neighbour-cell check exists. Skipping it to save nine lookups reintroduces exactly the boundary-miss failure the scheme is built to avoid.
• Precision that is too fine relative to drift causes cell-flapping, where the same physical position alternates between two or three adjacent cells fix to fix. If a downstream consumer treats a cell change as an event (entering/leaving a zone), flapping produces spurious enter/exit noise; debounce zone membership over a short window rather than trusting a single fix’s cell.
• The antimeridian and the poles are real edge cases, not theoretical ones, for any gateway fleet operating near ±180° longitude or in high latitudes. The wraparound in neighbour_cells handles longitude; latitude is clamped rather than wrapped because ±90° is a true boundary, not a seam — a zone that would need to wrap over a pole needs a different projection entirely, not a fix to this encoder.
• A geohash cell is not square except near the equator. Because latitude and longitude degrees do not correspond to equal ground distances outside the equator, the same bit precision yields a taller-than-wide or wider-than-tall cell at higher latitudes; pick precision from measured cell dimensions at the deployment’s actual latitude, not the equatorial table above taken at face value.
• Bucket hits are not containment. A fix landing in a cell registered to a zone is a candidate, not a confirmed match, whenever the zone’s true boundary does not align with cell edges — which is always, since cells are axis-aligned squares and zones rarely are. Treat lookup_zones results as the coarse pass that feeds a precise predicate, the same two-stage discipline used for the polygon containment check, whenever the zone boundary’s exact shape matters downstream.
• Coordinate precision upstream matters as much as bucketing precision downstream. Quantizing fixes to a fixed integer microdegree grid before encoding, as covered in spatial data precision standards, removes float-representation jitter that would otherwise occasionally flip a fix at a cell boundary for reasons unrelated to actual GPS drift.

Integrating with the reference dictionary #

The zone dictionary itself is built once at provisioning time by walking each zone’s bounding geometry and registering every cell its extent overlaps — structurally the same “register into every covered cell” step used to build the custom grid in the parent guide, just keyed by geohash integers instead of a bespoke (cx, cy) tuple. Persisting that dictionary as an integer-keyed store rather than rebuilding it every boot is exactly the staging question addressed in LMDB vs SQLite for spatial join staging; an LMDB environment keyed by the packed integer cell is a natural fit since the key is already a fixed-width machine word.

def register_zone(zone_dict: dict[int, list[str]], zone_id: str,
                   min_lon: float, min_lat: float, max_lon: float, max_lat: float,
                   bits: int = PRECISION_BITS):
    """Register a zone's bounding extent into every cell it overlaps, at load time."""
    _, _, step_lon, step_lat = cell_bounds(encode_cell(min_lat, min_lon, bits), bits)
    lon = min_lon
    while lon <= max_lon:
        lat = min_lat
        while lat <= max_lat:
            key = encode_cell(lat, lon, bits)
            zone_dict.setdefault(key, []).append(zone_id)
            lat += (step_lat - min_lat) if step_lat > min_lat else 0.01
        lon += (step_lon - min_lon) if step_lon > min_lon else 0.01

Related #

• Spatial Joins in Constrained Environments — the streaming grid-and-predicate join this bucketing scheme is an alternative front end for.
• LMDB vs SQLite for Spatial Join Staging — persisting the zone dictionary across restarts without a full rebuild.
• Spatial Data Precision Standards — quantizing fixes before encoding so boundary comparisons stay exact.