Skip to content

Geospatial indexing

TL;DR

  • An ordinary index sorts one dimension; "within 2 km" has two, so geospatial indexes flatten the map into cells or boxes a sorted index can scan.
  • Every proximity query is one pattern: cover the radius with cells, fetch candidates, filter by exact haversine distance, rank.
  • Geohash for key-value stores, an in-memory quadtree for skewed density, PostGIS for static relational data, H3 when neighbours must be equidistant.
  • Expect it in Yelp, Uber and Nearby Friends within five minutes.

Core concepts

Every scheme below does one thing: turn two coordinates into one sortable key, or into a tree whose boxes prune the search.

Latitude, longitude and haversine

A degree of latitude is ~111 km everywhere (40,000 km of circumference / 360); a degree of longitude shrinks with the cosine of the latitude: 111 km at the equator, 88 km in San Francisco, 56 km at 60 degrees north. Treating degrees as a flat plane therefore stretches east-west distances by 1/cos(lat), 27% at San Francisco, and breaks at the antimeridian, where 179.9 and -179.9 are 0.2 x 111 = 22 km apart, not 40,000. Haversine returns the great-circle distance in a handful of trigonometric calls (San Francisco to New York: 4,129 km) and is the exact filter every scheme below ends with.

Bounding-box queries

The index-friendly question is a rectangle: lat BETWEEN a AND b AND lon BETWEEN c AND d. A B-tree on latitude scans the whole band and filters longitude afterwards — 500k rows at 50k indexed reads/s is 10 s on one primary — and a composite index on (lat, lon) does not help: longitude is sorted only within equal latitudes. A circle of radius r fits the box lat +/- r/111 km, lon +/- r/(111 km x cos(lat)): every scheme below narrows to the box, then applies haversine.

Geohash: base32 cells, the precision table, neighbours and the boundary problem

A geohash interleaves bisections of longitude and latitude, longitude first: bit 1 is "east of 0 degrees", bit 2 "north of 0", bit 3 "east of the midpoint of what is left". Every five bits become one base32 character (digits and letters without a, i, l, o), so each extra character splits a cell 32 ways and a prefix is a cell: everything in 9q8yyk lies inside 9q8yy, the cell containing San Francisco.

Geohash grid

Derive the precision table rather than memorising it: precision p carries 5p bits, ceil(5p/2) for longitude and floor(5p/2) for latitude, so a cell is 40,000 km / 2^lon_bits wide and 20,000 km / 2^lat_bits tall.

Precision Cell at the equator (width x height) Resolves
1 5,004 km x 5,004 km a continent
2 1,251 km x 625 km a country
3 156 km x 156 km a region
4 39 km x 20 km a metro area
5 4.9 km x 4.9 km a district
6 1.2 km x 611 m a neighbourhood
7 153 m x 153 m a street
8 38 m x 19 m a building
9 4.8 m x 4.8 m a parking space

Odd precisions give square cells, even ones are twice as wide as tall, and width shrinks by cos(latitude): precision 5 is 3.9 km wide in San Francisco. Precision 6 tiles the planet with 32^6 = 2^30 cells — a cheap key in any store.

Two properties carry the design. Prefix sharing makes "the points in this cell" a prefix scan on any sorted index (geohash LIKE '9q8yy%'). The boundary problem is the catch: two points 22 m apart across the 45th parallel encode as c20bh2 and 9rbzur, sharing no character, because a bisection line runs between them. A shared prefix implies nearness, not the reverse, so every proximity query reads the query cell and its eight neighbours. Neighbours are arithmetic — cells at one precision are equal in degrees, so step the centre by one width or height and re-encode, wrapping at the antimeridian.

Quadtree

A quadtree is the adaptive alternative: a square that recursively splits into four quadrants whenever an insert would push a leaf past its capacity. Dense areas get deep and empty ones stay shallow, so a city and the ocean beside it cost what they contain, not a grid's worth of empty cells. A range query descends only into nodes whose box intersects the query rectangle: in the figure, 60 points at capacity 4 build 49 nodes and a query box visits 10 of them for its 6 points.

Quadtree

Nearest-neighbour search reuses the boxes with a priority queue: pop the closest item, expanding nodes into children keyed by the distance to their box and leaves into points keyed by exact distance; the k-th point popped is the k-th nearest, because everything still queued is farther. Sizing is easy: 200M businesses x 24 B (two 8 B coordinates and an 8 B id) = 4.8 GB, inside one server's 64-512 GB, so build it in memory at startup and replicate it as a read-only service. Its costs: the rebuild, a lock or copy-on-write for updates, and no stable cell id to store or shard by — which is why moving objects get a grid and static ones a tree.

R-trees

An R-tree is the balanced, disk-friendly tree of bounding rectangles behind PostGIS's GiST indexes and the spatial indexes of MySQL and SQLite: each node holds the smallest rectangle enclosing its children and an insert picks the subtree whose rectangle grows least. Rectangles overlap, so a search may follow several paths — the price of being the only family here that indexes shapes (roads, delivery zones, polygons) rather than points.

S2 and H3

Both fix a lat/lon grid's flaws. S2 (Google) projects the sphere onto the six faces of a cube, subdivides each face as a quadtree for 30 levels and orders the leaves along a Hilbert curve, so a 64-bit cell id's numeric neighbours are usually spatial neighbours, unlike geohash's Z-order, and cells at one level are near-uniform with no polar distortion. Its region coverer expresses any shape as a few mixed-level cells, turning "inside this polygon" into id range scans (Google Maps, MongoDB's 2dsphere). H3 (Uber) tiles the globe with hexagons at 16 resolutions (12 pentagons are unavoidable): every neighbour is equidistant, so a k-ring is a near-circle, where a square grid's diagonals sit 41% farther than its edges. Its hierarchy is approximate — a child is not exactly inside its parent — so it is an aggregation grid, not a storage key.

Redis GEO and PostGIS

Redis GEO stores members in a sorted set scored by a 52-bit geohash integer: GEOADD drivers lon lat id is one sorted-set insert and GEOSEARCH ... BYRADIUS 2 km runs the whole pattern below, inside one ~100k ops/s instance. 100k online drivers reporting every 4 s is 100k / 4 = 25k writes/s, a quarter of that — but one key lives on one shard, so key the sets by city and query the regions the radius touches. PostGIS adds a geography column with a GiST (R-tree) index: ST_DWithin(location, point, 2000) uses the index for the box and an exact check for the circle, and <-> returns k nearest neighbours. One primary serves 50k+ indexed reads/s before replicas: the right home for businesses, listings and polygons, the wrong one for 25k location updates a second.

The proximity query pattern

Whatever the index, the pattern is the flowchart below: the candidate set is always a few cells, never the table. Two things it cannot show: a fixed-precision index widens rather than choosing a precision (rings = ceil(r / cell size), 81 cells for 2 km at precision 6), and storing each moving object's current cell id makes a move "remove from the old cell, add to the new" — free when it stays put.

Cover the radius with cells, fetch candidates, filter exactly, then rank.

flowchart LR
    q["Query: lat, lon, radius r, limit K"] --> prec["Precision p: the finest cell<br/>still at least r wide and tall"]
    prec --> cover["Covering cells: query cell<br/>plus its 8 neighbours"]
    cover --> idx[("Per-cell lookup<br/>Redis GEO, cell to ids map,<br/>or a prefix scan on geohash")]
    idx --> cand["Candidate ids<br/>(cells are rectangles, the circle is not)"]
    cand --> exact["Exact filter:<br/>haversine(point, query) within r"]
    exact --> rank["Sort by distance, cut at K,<br/>join details"]
    rank --> more{"Fewer than K hits?"}
    more -->|"yes"| widen["Add one ring of cells"]
    widen --> idx
    more -->|"no"| resp["Nearest K"]

Pick the index from what moves and what is queried.

flowchart TD
    q_static{"Static objects already in a<br/>relational database (shops, listings)?"} -->|"yes"| pg["PostGIS: GiST R-tree,<br/>ST_DWithin and the KNN operator"]
    q_static -->|"no"| q_move{"Objects move constantly<br/>(drivers, couriers, friends)?"}
    q_move -->|"yes"| q_hex{"Neighbours must be equidistant<br/>(surge zones, smoothing)?"}
    q_hex -->|"yes"| h3["H3: a hexagon id per update,<br/>k-ring queries"]
    q_hex -->|"no"| redis["Redis GEO or a geohash map<br/>sharded by region"]
    q_move -->|"no"| q_skew{"Density varies by orders<br/>of magnitude (cities vs oceans)?"}
    q_skew -->|"yes"| quad["In-memory quadtree<br/>or an S2 covering"]
    q_skew -->|"no"| geo["Geohash prefix index<br/>in any sorted store"]
    pg --> filt["Always: exact haversine filter, then rank"]
    h3 --> filt
    redis --> filt
    quad --> filt
    geo --> filt

Trade-offs

Index Cells Adapts to density Neighbours Moving objects Shapes Typical home
Geohash Rectangles, 2:1 at even precisions No 8 by arithmetic; boundary problem Cheap: re-encode, move one key No Key-value stores, Redis GEO, prefix indexes
Quadtree Squares split on demand Yes Tree walk Delete and reinsert under a lock Points only In-memory service, rebuilt from a snapshot
R-tree Bounding rectangles Yes Overlapping paths Moderate Yes PostGIS, MySQL, SQLite
S2 Near-uniform squares on a cube, Hilbert-ordered Mixed-level coverings Cell id arithmetic Cheap Via coverings Google Maps, MongoDB 2dsphere
H3 Hexagons (12 pentagons) 16 resolutions Equidistant k-ring Cheap Approximate Uber analytics, surge, ETA buckets

Choose from the write rate and the shape of the data. Static points already in a relational database belong in PostGIS: one index, exact results, polygons when you need them. Objects that move every few seconds want a cell key you can overwrite cheaply — a geohash or a Redis GEO set keyed by region — plus nine cells to defuse the boundary problem. When density is wildly uneven and the set fits in memory, a quadtree answers range and nearest-K queries in a handful of node visits with no empty cells, paid for with a rebuild on restart and a lock on updates: snapshot it periodically, treat it as read-mostly. Reach for S2 for polygons over a spherical grid at global scale, and for H3 when the question is about neighbourhoods rather than points. Whatever you pick, the exact distance filter is not optional: the index narrows, it never decides.

Python implementation

encode bisects longitude and latitude alternately, five bits per base32 character; bounds replays those bits and decode returns the cell centre:

code/hld/geohash.py — encode, bounds, decode
@dataclass(frozen=True, slots=True)
class Bounds:
    """The rectangle a geohash names, in degrees."""

    lat_lo: float
    lat_hi: float
    lon_lo: float
    lon_hi: float

    @property
    def center(self) -> tuple[float, float]:
        return (self.lat_lo + self.lat_hi) / 2, (self.lon_lo + self.lon_hi) / 2

    @property
    def height_deg(self) -> float:
        return self.lat_hi - self.lat_lo

    @property
    def width_deg(self) -> float:
        return self.lon_hi - self.lon_lo

    def contains(self, lat: float, lon: float) -> bool:
        return self.lat_lo <= lat < self.lat_hi and self.lon_lo <= lon < self.lon_hi


def encode(lat: float, lon: float, precision: int = 9) -> str:
    """Geohash of a point: bisect longitude, then latitude, then longitude... and emit five bits
    per base32 character. ``precision`` characters carry 5 x precision bits."""
    _validate_point(lat, lon)
    if not 1 <= precision <= MAX_PRECISION:
        raise ValidationError(f"precision must be 1-{MAX_PRECISION}")
    lat_lo, lat_hi, lon_lo, lon_hi = -90.0, 90.0, -180.0, 180.0
    chars: list[str] = []
    bits = bit_count = 0
    longitude_turn = True
    while len(chars) < precision:
        if longitude_turn:
            mid = (lon_lo + lon_hi) / 2
            if lon >= mid:
                bits, lon_lo = bits * 2 + 1, mid
            else:
                bits, lon_hi = bits * 2, mid
        else:
            mid = (lat_lo + lat_hi) / 2
            if lat >= mid:
                bits, lat_lo = bits * 2 + 1, mid
            else:
                bits, lat_hi = bits * 2, mid
        longitude_turn = not longitude_turn
        bit_count += 1
        if bit_count == 5:
            chars.append(BASE32[bits])
            bits = bit_count = 0
    return "".join(chars)


def bounds(geohash: str) -> Bounds:
    """Replay the bisections a geohash encodes to recover its cell."""
    _validate_hash(geohash)
    lat_lo, lat_hi, lon_lo, lon_hi = -90.0, 90.0, -180.0, 180.0
    longitude_turn = True
    for ch in geohash:
        value = _DECODE[ch]
        for shift in (4, 3, 2, 1, 0):
            bit = (value >> shift) & 1
            if longitude_turn:
                mid = (lon_lo + lon_hi) / 2
                lon_lo, lon_hi = (mid, lon_hi) if bit else (lon_lo, mid)
            else:
                mid = (lat_lo + lat_hi) / 2
                lat_lo, lat_hi = (mid, lat_hi) if bit else (lat_lo, mid)
            longitude_turn = not longitude_turn
    return Bounds(lat_lo, lat_hi, lon_lo, lon_hi)


def decode(geohash: str) -> tuple[float, float]:
    """Centre of the cell: the best point estimate, off by at most half a cell in each axis."""
    return bounds(geohash).center

adjacent steps one cell in any of eight directions; cell_size_km, precision_for_radius_km and cells_covering turn radii into precisions and coverings:

code/hld/geohash.py — neighbours, cell sizes, coverings
class Direction(Enum):
    """Steps in (latitude, longitude) cell units, clockwise from north."""

    N = (1, 0)
    NE = (1, 1)
    E = (0, 1)
    SE = (-1, 1)
    S = (-1, 0)
    SW = (-1, -1)
    W = (0, -1)
    NW = (1, -1)


def adjacent(geohash: str, direction: Direction) -> str | None:
    """The cell one step away, or ``None`` beyond a pole.

    Every cell at a given precision has the same size in degrees, so stepping from the centre
    by one cell height or width and re-encoding is exact and sidesteps the base32 irregularity.
    Longitude wraps at the antimeridian.
    """
    box = bounds(geohash)
    d_lat, d_lon = direction.value
    lat = box.center[0] + d_lat * box.height_deg
    lon = box.center[1] + d_lon * box.width_deg
    if not -90 <= lat <= 90:
        return None
    if lon >= 180:
        lon -= 360
    elif lon < -180:
        lon += 360
    return encode(lat, lon, len(geohash))


def neighbors(geohash: str) -> list[str]:
    """The surrounding cells clockwise from north: eight, fewer at a pole."""
    return [cell for d in Direction if (cell := adjacent(geohash, d)) is not None]


def cell_size_km(precision: int, lat: float = 0.0) -> tuple[float, float]:
    """``(width, height)`` of a cell in km at latitude ``lat``.

    5 x precision bits split longitude-first, so longitude gets ceil(5p/2) bits and latitude
    floor(5p/2): odd precisions give square cells, even ones are twice as wide as tall. Width
    shrinks by cos(lat) away from the equator; height does not.
    """
    if not 1 <= precision <= MAX_PRECISION:
        raise ValidationError(f"precision must be 1-{MAX_PRECISION}")
    lon_bits = (5 * precision + 1) // 2
    lat_bits = 5 * precision // 2
    width = 2 * math.pi * EARTH_RADIUS_KM / 2**lon_bits * math.cos(math.radians(lat))
    height = math.pi * EARTH_RADIUS_KM / 2**lat_bits
    return width, height


def precision_for_radius_km(radius_km: float, lat: float = 0.0) -> int:
    """The finest precision whose cell is still at least ``radius_km`` in both directions, so
    the query cell plus its eight neighbours contain every point within the radius."""
    if radius_km <= 0:
        raise ValidationError("radius must be positive")
    for precision in range(MAX_PRECISION, 1, -1):
        if min(cell_size_km(precision, lat)) >= radius_km:
            return precision
    return 1


def cells_covering(lat: float, lon: float, radius_km: float, precision: int) -> list[str]:
    """Cells at ``precision`` whose union contains every point within ``radius_km``: the query
    cell plus ``rings`` rings around it, where one ring is the classic nine-cell search."""
    _validate_point(lat, lon)
    if radius_km <= 0:
        raise ValidationError("radius must be positive")
    far_lat = min(90.0, abs(lat) + radius_km / KM_PER_DEGREE)  # narrowest cells the circle touches
    rings = max(1, math.ceil(radius_km / min(cell_size_km(precision, far_lat))))
    if rings > MAX_RINGS:
        raise ValidationError(
            f"a {radius_km} km radius needs {rings} rings of precision-{precision} cells; "
            "index at a coarser precision"
        )
    box = bounds(encode(lat, lon, precision))
    c_lat, c_lon = box.center
    cells: list[str] = []
    for d_lat in range(-rings, rings + 1):
        step_lat = c_lat + d_lat * box.height_deg
        if not -90 <= step_lat <= 90:
            continue
        for d_lon in range(-rings, rings + 1):
            step_lon = (c_lon + d_lon * box.width_deg + 180) % 360 - 180
            cell = encode(step_lat, step_lon, precision)
            if cell not in cells:
                cells.append(cell)
    return cells

haversine_km is the exact filter:

code/hld/geohash.py — haversine
def haversine_km(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    """Great-circle distance. Euclidean distance on degrees is wrong twice over: a degree of
    longitude shrinks with cos(latitude), and the surface curves."""
    _validate_point(lat1, lon1)
    _validate_point(lat2, lon2)
    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    d_phi = phi2 - phi1
    d_lambda = math.radians(lon2 - lon1)
    a = math.sin(d_phi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(d_lambda / 2) ** 2
    return 2 * EARTH_RADIUS_KM * math.asin(math.sqrt(a))

GeoIndex is the query pattern itself: cell-to-ids and id-to-point maps under one lock, then candidates, haversine, sort, limit:

code/hld/geohash.py — the proximity index
@dataclass(frozen=True, slots=True)
class Hit:
    item_id: str
    lat: float
    lon: float
    distance_km: float


class GeoIndex:
    """Proximity search over geohash cells at one fixed precision.

    ``_cells`` maps a cell to the ids inside it and ``_points`` maps an id to its coordinates and
    cell; ``_lock`` guards both, because items move (drivers, couriers) while queries run. The
    query pattern: covering cells, candidate ids, exact haversine filter, sort by distance.
    """

    def __init__(self, precision: int = 6) -> None:
        if not 1 <= precision <= MAX_PRECISION:
            raise ValidationError(f"precision must be 1-{MAX_PRECISION}")
        self._precision = precision
        self._cells: dict[str, set[str]] = {}
        self._points: dict[str, tuple[float, float, str]] = {}
        self._lock = threading.Lock()

    def __len__(self) -> int:
        with self._lock:
            return len(self._points)

    @property
    def precision(self) -> int:
        return self._precision

    def add(self, item_id: str, lat: float, lon: float) -> str:
        """Insert an item or move it; returns the cell it now lives in."""
        if not item_id:
            raise ValidationError("item_id must be non-empty")
        _validate_point(lat, lon)
        cell = encode(lat, lon, self._precision)
        with self._lock:
            self._unlink(item_id)
            self._points[item_id] = (lat, lon, cell)
            self._cells.setdefault(cell, set()).add(item_id)
        return cell

    def remove(self, item_id: str) -> None:
        with self._lock:
            if item_id not in self._points:
                raise NotFoundError(f"unknown item {item_id!r}")
            self._unlink(item_id)

    def _unlink(self, item_id: str) -> None:
        """Detach an item from its cell; caller holds the lock."""
        previous = self._points.pop(item_id, None)
        if previous is None:
            return
        members = self._cells[previous[2]]
        members.discard(item_id)
        if not members:
            del self._cells[previous[2]]

    def nearby(self, lat: float, lon: float, radius_km: float, limit: int = 10) -> list[Hit]:
        """Items within ``radius_km``, nearest first, at most ``limit`` of them."""
        if limit <= 0:
            raise ValidationError("limit must be positive")
        cells = cells_covering(lat, lon, radius_km, self._precision)
        with self._lock:
            candidates = [
                (item_id, *self._points[item_id][:2])
                for cell in cells
                for item_id in self._cells.get(cell, ())
            ]
        hits = [
            Hit(item_id, p_lat, p_lon, distance)
            for item_id, p_lat, p_lon in candidates
            if (distance := haversine_km(lat, lon, p_lat, p_lon)) <= radius_km
        ]
        hits.sort(key=lambda hit: (hit.distance_km, hit.item_id))
        return hits[:limit]

    def candidates(self, lat: float, lon: float, radius_km: float) -> tuple[list[str], int]:
        """``(covering cells, candidate count)`` for a query: the cost before the exact filter."""
        cells = cells_covering(lat, lon, radius_km, self._precision)
        with self._lock:
            count = sum(len(self._cells.get(cell, ())) for cell in cells)
        return cells, count

uv run python -m hld.geohash prints:

encode(37.7749, -122.4194) by precision: 9 | 9q | 9q8 | 9q8y | 9q8yy | 9q8yyk | 9q8yyk8 | 9q8yyk8y
decode('9q8yy') = (37.7710, -122.4097); the cell is 3.9 km x 4.9 km, so the error is at most half of that
precision table, width x height at the equator (width shrinks by cos(lat): x0.79 at SF):
  1: 5,004 km x 5,004 km     5: 4.9 km x 4.9 km
  2: 1,251 km x 625 km       6: 1.2 km x 611 m
  3: 156 km x 156 km         7: 153 m x 153 m
  4: 39 km x 20 km           8: 38.2 m x 19.1 m
neighbours of 9q8yy: N=9q8zn, NE=9q8zp, E=9q8yz, SE=9q8yx, S=9q8yw, SW=9q8yt, W=9q8yv, NW=9q8zj
boundary problem: c20bh2 and 9rbzur are 22 m apart and share no prefix
search 2 km around Union Square: precision 5 (3.9 km x 4.9 km cells), 9 cells, 3 of 7 places are candidates
   0.4 km  Powell St station
   1.5 km  Ferry Building
   1.6 km  Coit Tower
  never read (outside the 9 cells): Golden Gate Bridge (7 km), Oakland, Jack London Sq (11 km), Berkeley campus (16 km), Palo Alto (45 km)

Rect distinguishes half-open node boxes from closed query boxes and measures point-to-box distance, which nearest orders its heap by:

code/hld/quadtree.py — points and boxes
@dataclass(frozen=True, slots=True)
class Point:
    x: float
    y: float
    item_id: str = ""


@dataclass(frozen=True, slots=True)
class Rect:
    """Axis-aligned box. Node boxes are half-open (a point on the shared edge belongs to one
    child only); query boxes are closed, so a point on the query edge counts."""

    x0: float
    y0: float
    x1: float
    y1: float

    def __post_init__(self) -> None:
        if self.x1 <= self.x0 or self.y1 <= self.y0:
            raise ValidationError("a Rect needs x0 < x1 and y0 < y1")

    @property
    def width(self) -> float:
        return self.x1 - self.x0

    @property
    def height(self) -> float:
        return self.y1 - self.y0

    def contains(self, p: Point) -> bool:
        return self.x0 <= p.x < self.x1 and self.y0 <= p.y < self.y1

    def covers(self, p: Point) -> bool:
        return self.x0 <= p.x <= self.x1 and self.y0 <= p.y <= self.y1

    def intersects(self, other: Rect) -> bool:
        return not (
            self.x0 > other.x1 or self.x1 < other.x0 or self.y0 > other.y1 or self.y1 < other.y0
        )

    def distance_sq(self, x: float, y: float) -> float:
        """Squared distance from a point to the nearest point of the box (0 inside it)."""
        dx = max(self.x0 - x, 0.0, x - self.x1)
        dy = max(self.y0 - y, 0.0, y - self.y1)
        return dx * dx + dy * dy

QuadTree.insert splits a full leaf on the way down, query prunes by box intersection and counts nodes visited, and nearest is a best-first search over one heap of nodes and points:

code/hld/quadtree.py — insert, range query, nearest K
@dataclass(frozen=True, slots=True)
class RangeResult:
    points: list[Point]
    nodes_visited: int


class QuadTree:
    """Points in a square region, split adaptively; ``_lock`` guards the whole tree."""

    def __init__(self, boundary: Rect, capacity: int = 4, max_depth: int = 16) -> None:
        if capacity <= 0 or max_depth <= 0:
            raise ValidationError("capacity and max_depth must be positive")
        self._root = _Node(boundary, depth=0)
        self._capacity = capacity
        self._max_depth = max_depth
        self._size = 0
        self._lock = threading.Lock()

    def __len__(self) -> int:
        with self._lock:
            return self._size

    @property
    def boundary(self) -> Rect:
        return self._root.box

    def insert(self, p: Point) -> None:
        """Walk down to the leaf for ``p``; split it while it is full, unless ``max_depth`` is
        reached (which is what stops duplicate coordinates from splitting forever)."""
        if not self._root.box.contains(p):
            raise ValidationError(f"{p} is outside the tree boundary {self._root.box}")
        with self._lock:
            node = self._root
            while node.children is not None:
                node = node.child_for(p)
            while len(node.points) >= self._capacity and node.depth < self._max_depth:
                node.subdivide()
                node = node.child_for(p)
            node.points.append(p)
            self._size += 1

    def query(self, rect: Rect) -> RangeResult:
        """Every point inside ``rect``; only nodes whose box meets ``rect`` are visited."""
        found: list[Point] = []
        visited = 0
        with self._lock:
            stack = [self._root]
            while stack:
                node = stack.pop()
                if not node.box.intersects(rect):
                    continue
                visited += 1
                if node.children is None:
                    found.extend(p for p in node.points if rect.covers(p))
                else:
                    stack.extend(node.children)
        return RangeResult(found, visited)

    def nearest(self, x: float, y: float, k: int = 1) -> list[tuple[float, Point]]:
        """The ``k`` closest points as ``(distance, point)``, nearest first.

        One heap holds nodes (keyed by distance to their box) and points (exact distance). A
        point popped before any node is closer than anything those nodes could contain.
        """
        if k <= 0:
            raise ValidationError("k must be positive")
        tiebreak = itertools.count()
        heap: list[tuple[float, int, _Node | Point]] = []
        result: list[tuple[float, Point]] = []
        with self._lock:
            heapq.heappush(heap, (self._root.box.distance_sq(x, y), next(tiebreak), self._root))
            while heap and len(result) < k:
                dist_sq, _, item = heapq.heappop(heap)
                if isinstance(item, Point):
                    result.append((dist_sq**0.5, item))
                elif item.children is None:
                    for p in item.points:
                        d = (p.x - x) ** 2 + (p.y - y) ** 2
                        heapq.heappush(heap, (d, next(tiebreak), p))
                else:
                    for child in item.children:
                        heapq.heappush(heap, (child.box.distance_sq(x, y), next(tiebreak), child))
        return result

    def stats(self) -> tuple[int, int]:
        """``(node count, max depth)``."""
        with self._lock:
            count = depth = 0
            stack = [self._root]
            while stack:
                node = stack.pop()
                count += 1
                depth = max(depth, node.depth)
                stack.extend(node.children or [])
        return count, depth

uv run python -m hld.quadtree prints:

60 seeded points in a 100 x 100 box, capacity 4: 49 nodes, depth 3
range query x 60-90, y 15-45: 6 points, 10 of 49 nodes visited (a scan would test all 60 points)
  (60.9, 15.3), (63.6, 36.5), (65.5, 39.6), (79.2, 42.2), (87.6, 26.3), (87.6, 31.5)
nearest 3 to (50, 50):
  p17  (37.9, 55.2) at 13.2
  p7   (65.0, 54.5) at 15.6
  p27  (64.8, 60.9) at 18.4
200 more points inside a 2 x 2 patch around (10, 10): 217 nodes, depth 10; only that branch got deeper
same range query elsewhere: 6 points, 10 nodes visited; the patch itself: 201 points, 151 nodes visited
nearest to (10, 10): c10 (dense branch, still a heap walk)

In the interview

Say where coordinates live and how they are queried when you reach the data model: "Drivers publish a location every 4 s into a Redis GEO set keyed by city. Matching asks for drivers within 2 km: the query cell plus its eight neighbours, an exact distance filter, the nearest K." That covers index, sharding, the boundary problem and ranking in one sentence.

Phrases that signal depth: "the cell plus its eight neighbours, because of the boundary problem"; "precision from the radius"; "the index returns candidates, haversine returns answers".

Why not index latitude and longitude as two ordinary columns?

A B-tree sorts one key. An index on (lat, lon) scans the whole latitude band and filters longitude afterwards: 10 s for 500k rows at 50k indexed reads/s. Cells and boxes prune both dimensions at once.

A driver sits on a cell boundary. Do you miss them?

No: the covering is the query cell plus its eight neighbours, arithmetic on the cell centre, so the extra cost is eight key lookups and haversine drops what falls outside the radius.

The first ring returns three drivers but the rider wants ten. What now?

Widen by one ring at the same precision (25 cells, then 49) rather than jumping to the city; the exact filter keeps results correct. A quadtree does nearest-K natively.

How do you shard a global proximity index?

By region: a geohash prefix or city id is the partition key, so a query touches one shard unless the radius crosses a border and fans out to two. Hot cities get their own shard, with finer cells.

Drivers move every few seconds. How do you keep the index fresh without thrashing it?

GEOADD is an upsert, one write per update: 100k drivers every 4 s is 25k writes/s. Store the driver's current cell so an update inside it is a no-op, and set a TTL so a silent driver disappears.

Interview tip

Derive one cell size aloud instead of quoting a table: "precision 6 is 15 longitude bits, 40,000 km / 2^15 = 1.2 km wide". It shows you understand the encoding, not the table.

Common mistakes

  • Reading only the query cell: results vanish for users near a cell edge. Fix: cell plus eight neighbours, then the exact filter.
  • Euclidean distance on degrees: east-west distances inflate by 1/cos(lat), 27% in San Francisco, and the antimeridian becomes a 40,000 km gap. Fix: haversine, planar only inside a city-sized window.
  • A precision chosen without the radius: precision 8 needs thousands of 38 m x 19 m cells for a 2 km circle; precision 3 reads a 156 km region. Fix: the finest precision whose cell is at least r.
  • The wrong home for moving objects: 25k updates/s against a GiST index on one primary, or every driver on the planet in one Redis sorted set capped at ~100k ops/s. Fix: a cell key in Redis or an in-memory grid keyed by region, with the relational store for static attributes.

Common mistake

Returning the cells' contents as the answer. Cells are rectangles and the circle is not: nine cells around a 2 km query at precision 5 cover 9 x 3.9 x 4.9 = 172 km², fourteen times the 12.6 km² circle, so most of what you return is outside the radius. The index narrows; the distance decides.

Self-check

How large is a precision-6 geohash cell, and how many are there?

15 bits each way: 40,000 km / 2^15 = 1.2 km wide and 20,000 km / 2^15 = 611 m tall at the equator; 32^6 = 2^30 cells.

Why are even-precision cells twice as wide as they are tall?

Bits alternate longitude first, so an even precision splits both axes equally while longitude spans 360 degrees against latitude's 180. An odd precision gives longitude one extra split, squaring the cell.

What is the pruning rule in a quadtree range query?

Descend into a node only if its box intersects the query rectangle; at a leaf, test its points. Everything outside is skipped.

Why does Uber use hexagons for surge pricing?

Every neighbour of a hexagon is equidistant, so a k-ring is a near-circle and smoothing treats neighbours equally; square-grid diagonals are 41% farther.

What does Redis store when you call GEOADD?

A sorted-set member scored by a 52-bit interleaved geohash; a radius search scans the covering cells' score ranges and filters by exact distance.