Geospatial and Time-Series Data Questions

Direct answer

For city-scale real-time matching the workload is write-heavy: every online driver sends a location every few seconds, while ride requests are far rarer. That favours structures with O(1) updates: a uniform grid or geohash grid (cell ID to a set of drivers) held in memory and sharded by city or region. Quadtrees adapt to uneven density but pay tree maintenance on every move; KD-trees are excellent for static nearest-neighbour search but degrade under constant updates, so they fit only as periodically rebuilt snapshots; R-trees are built for rectangles and polygons on disk, so they belong to the static layer (airports, pickup zones, city boundaries), not to the moving drivers. My pick: an in-memory geohash (or H3) cell grid per city shard for live drivers, plus an R-tree for static geofences.

Terms used

• Uniform grid (grid-based tiling): the map is cut into equal cells; each cell keeps a list of the drivers inside it.
• Geohash grid: the same idea, with cells named by geohash strings so they nest (a cell's children share its prefix) and can be stored in a sorted key-value store.
• Quadtree: a tree where each node covers a square and splits into four children when it holds too many points, so dense areas get small cells and empty areas stay large.
• KD-tree: a binary tree that alternately splits points by x and by y at the median.
• R-tree: a tree of nested bounding rectangles, where each node's rectangle encloses its children; designed for shapes and for disk pages.
• k: the number of results returned by a query.
• Redis: an in-memory key-value data store. Its GEOADD/GEOSEARCH commands keep locations as geohash scores in a sorted set, a data structure that keeps its members ordered by a score so a range of scores can be fetched quickly.
• N and M (in the Redis complexity figures below): Redis's own docs define these per command, and they are not the same variable each time. For GEOADD, N is the number of elements already in the sorted set (the whole index). For GEOSEARCH, N is the number of elements in the grid cells that coarsely overlap the query shape, and M is the number of those elements that actually fall inside the exact shape once filtered. Both differ from this answer's k.

Comparison

Framing by workload characteristics

Worked numbers for one large city (example inputs): 50,000 online drivers pinging every 4 s, peak 300 ride requests per second.

• Update frequency: 50,000 / 4 = 12,500 location writes per second. Only those that cross a cell boundary change the index; with 150 m cells and cars averaging 8 m/s, a ping moves about 8 x 4 = 32 m, so on a simple model where the chance of landing in a new cell scales with (step length) / (cell size), that is 32 / 150, roughly 1 in 5: most pings stay in the same cell, and only about 12,500 x (32 / 150), around 2,700, writes per second actually touch the index.
• Query rate (QPS, queries per second): 300 ride requests/s. If each triggers a nearby-driver search plus a few retries (expanding the search radius), call it about 1,000 searches/s.
• Read/write ratio: roughly 1,000 : 12,500, about 1 : 12.5. Writes dominate, so the structure must make writes nearly free, which rules out anything that rebalances on insert.
• Data size: 50,000 drivers x (8-byte ID + two 8-byte coordinates + cell ID) is a few megabytes, so the whole index fits in one process's memory; durability is unnecessary because positions are rebuilt from the next ping within seconds.

Placement in a distributed real-time matching architecture

• Partition by geography: one shard (or a primary node that takes writes, plus a replica that mirrors it for redundancy and read scaling) per city or per region of a large city. Queries are local, so almost all traffic stays inside one shard.
• Boundary queries: a rider near a shard border needs cells owned by the neighbour shard; the matcher fans out (sends the same query to several shards in parallel) to at most the few adjacent shards and merges candidates.
• Hot spots: an airport or stadium cell can hold thousands of drivers. The grid handles this by switching those areas to a finer precision, or by storing drivers in the hot cell in a small secondary structure. A quadtree handles it automatically, which is its main advantage.
• Snapshots for heavy reads: if a batch matcher needs many k-nearest-neighbour queries at once, build a KD-tree snapshot of the grid every second or two and query that, accepting a slightly stale view and re-checking the chosen drivers' latest positions.
• Static layer: airports, restricted zones and pricing areas are polygons that change rarely: an R-tree (for example the GiST, Generalized Search Tree, index in PostGIS, the spatial extension for PostgreSQL) answers "which zone is this point in" quickly.

Recommendation and what would flip it

City-scale real-time matching: in-memory geohash (or H3) grid, sharded by city, cell size chosen so a typical search touches 9 to 25 cells (a 3x3 block, one ring of neighbours, at the tightest; a 5x5 block, two rings, when the first ring comes up short), with an exact distance or travel-time filter afterwards. I would flip to a quadtree if density varied by orders of magnitude within one shard and fixed cells produced either huge hot cells or many empty lookups. I would use a KD-tree only for offline or batch nearest-neighbour work over a static snapshot, and an R-tree when the indexed objects are shapes or the index must live on disk in a database.

Pitfalls

• Benchmarking with query-only workloads, which makes KD-trees look best while hiding the update cost.
• Choosing structures by big-O alone: at 50,000 points, constant factors, cache locality (how often the data a CPU needs is already in its fast on-chip cache rather than main memory) and lock contention (threads waiting on each other for the same lock) on hot cells matter more.
• Persisting every ping to a database index (for example an R-tree) at 12,500 writes/s per city, which turns an in-memory problem into a write-amplification problem (one logical write turning into many more physical disk writes than expected).

Direct answer

Store raw readings in a table partitioned by day and keep 30 days of partitions; build hourly rollups from raw continuously and keep them 3 years; build daily rollups from the hourly table and keep them 3 years too, purely to make multi-year queries fast. Expire data by dropping whole partitions, never by DELETE. Every query carries a time predicate so the database reads only the partitions in range (partition pruning). A thin query layer (a view or a small service) routes each request: the part of the range within the last 30 days comes from raw, the older part from the rollup tier, stitched at a fixed boundary. At 10 million points per day this keeps the database around 16 GB instead of more than 500 GB for 3 years of raw data.

Key terms

• Retention: how long a tier keeps data. SLA (service-level agreement) here means the commitment: raw at full resolution for 30 days, summaries for 3 years.
• Rollup: a pre-computed summary per sensor per bucket (hour, day): min, max, sum, count.
• Compaction: replacing many small units of data with fewer, larger ones. Here it means turning raw rows into rollup rows and dropping the raw rows, and, for columnar storage, merging many small files into large ones.
• Partition: a physical slice of a table, here one per day for raw data. The database treats it as a separate table underneath, so dropping it is instant.
• Partition pruning: the planner (the part of Postgres that decides how to execute a query, including which partitions it actually needs to open) skips partitions whose time range cannot match the WHERE clause. A query for "yesterday" opens 1 partition out of 30.

The tiering table

Storage cost vs latency, the explicit trade-off. Each coarser tier costs storage and one more job, and buys query latency on long ranges. The hourly tier is mandatory (the 3-year SLA requires summaries to exist once raw is gone). The daily tier is optional: it costs about 4% of the hourly tier's storage (1/24) and cuts a 3-year chart from 26,280 rows per sensor to 1,095. That is cheap latency, so I keep it. A further tier (weekly, monthly) would save almost nothing more and is not worth the extra job.

Worked example at 10 million points per day

Assume 1,000 sensors, each sending one reading about every 8.64 seconds (86,400 / 10,000), so 10,000,000 points per day. Take 48 bytes per raw row in a row store (timestamp, sensor id, value, plus per-row overhead) and 64 bytes per rollup row; these are planning assumptions, not measurements.

Total about 16.2 GB against 525.6 GB. Columnar compression would shrink both sides, but the ratio is what matters.

Implementing the rollups

• Hourly job, every 15 minutes, recomputes the last 3 hours of buckets from raw with an upsert (INSERT ... ON CONFLICT (sensor_id, hour) DO UPDATE). Recomputing a trailing window picks up readings that arrived late, for example from a sensor that reconnected. Upserts make the job idempotent: running it twice gives the same result, so retries after a crash are safe.
• Daily job runs after midnight plus a margin (say 03:00) and aggregates the previous day's 24 hourly rows: min of mins, max of maxes, sum of sums, sum of counts. These combine exactly, which an average of averages would not.
• Watermark table: each job records the latest bucket it completed. The retention job reads it (next section).

Retention and compaction by partition

-- PostgreSQL declarative partitioning, one partition per day for raw
CREATE TABLE readings (
  sensor_id int NOT NULL,
  ts        timestamptz NOT NULL,
  value     double precision NOT NULL
) PARTITION BY RANGE (ts);

CREATE TABLE readings_2026_09_28 PARTITION OF readings
  FOR VALUES FROM ('2026-09-28') TO ('2026-09-29');
CREATE INDEX ON readings (sensor_id, ts);  -- propagates to every partition

A scheduler (pg_partman, a PostgreSQL extension that creates and drops partitions on a schedule, or a cron job) creates partitions a few days ahead and drops old ones:

1. Check the hourly watermark covers the whole partition's day, plus the 3-hour late-data window the hourly job recomputes (see "Implementing the rollups" above), so a drop never runs ahead of data the rollup job might still need.
2. ALTER TABLE readings DETACH PARTITION readings_2026_08_28; then DROP TABLE readings_2026_08_28;
3. If the watermark is behind, do not drop: alert instead. A 31st day of raw data costs 0.48 GB; a lost day of summaries is unrecoverable.

Dropping a partition is a metadata operation (it just removes an entry from the table's list of partitions and unlinks a file; it never opens or touches a row inside it). DELETE FROM readings WHERE ts < now() - interval '30 days' would touch 10 million rows a day: each deleted row still has to be recorded in the write-ahead log (WAL, the append-only file Postgres writes every change to before applying it, so it can recover after a crash), and each becomes a dead row (an old row version Postgres keeps around rather than erasing in place) that VACUUM (Postgres's background cleanup process) later has to scan and reclaim.

Pruning only works if queries let it. WHERE ts >= '2026-09-01' AND ts < '2026-09-08' prunes to 7 partitions. WHERE date_trunc('day', ts) = '2026-09-01' wraps the column in a function and may force a scan of all 30.

Serving queries that need both raw and rolled-up data

flowchart LR
  C[Client: sensor, range, resolution] --> R{Query router}
  R -->|range within 30 days| RAW[(Raw partitions)]
  R -->|older part of range| H[(Hourly)]
  R -->|range over 1 year| D[(Daily)]
  RAW --> M[Merge at boundary]
  H --> M
  D --> M
  M --> C

• Rule 1: choose the output resolution first. A 90-day chart at hourly resolution is 2,160 points per sensor. Serve it entirely from the hourly table (which covers the last 30 days too), rather than stitching, so the line is consistent.
• Rule 2: stitch only when the caller needs raw detail for the recent part. Example: "last 60 days, raw where available". The router splits at the boundary now() - 30 days rounded down to a whole hour: raw rows after it, hourly rows before it, returned with a resolution field on each point so the client can draw them differently.
• Rule 3: return the same shape from every tier (min, max, avg as sum / count, count), so a raw point is presented as a bucket of count 1.
• A SQL UNION ALL view across raw and hourly with complementary time predicates works for simple cases; a service layer is better once you have several tiers and caching.

Trade-offs and pitfalls

• Dropping raw before the rollup covers it is the one failure that loses data permanently. Gate every drop on the watermark.
• Percentiles cannot be rebuilt from min, max, sum, count. If the 3-year SLA includes percentiles, store a fixed-bucket histogram per hour (a count of readings falling into fixed value ranges, like 0-10, 10-20, which can be summed across hours and used to estimate a percentile, unlike a single stored percentile number).
• Too-small partitions (hourly partitions for raw) mean 720 partitions to plan over for 30 days and slower planning; too-large ones (monthly) mean retention can only be exact to the month. Daily matches the 30-day SLA.
• What would change it: at 1 billion points per day instead of 10 million, raw storage moves to a columnar TSDB (a time-series database that stores each column of data compressed together, built for that scale) or Parquet files (a compressed, columnar file format) in object storage, but the tier table, watermark-gated drops and the router stay the same.

Direct answer

Store (or index) each listing as a geography point, put a partial GiST index on it covering only active listings, and drive the per-search lookup with a LEFT JOIN LATERAL subquery that filters with ST_DWithin(..., 5000) and orders by the <-> distance operator with LIMIT 5. The lateral join runs one small index probe per search instead of comparing every search against every listing. The city percentage is a second aggregation over the per-search result counts. One gap in the stated schema has to be closed first: searches has no city column and no key, so the answer adds a search_id and assigns a city with a point-in-polygon join (a join whose condition is "this point falls inside this polygon", here testing whether a search's coordinate lands inside a city's boundary shape) against a cities boundary table.

Why each piece is there

Key terms.

• A spatial index lets the database skip rows that are nowhere near the query point. PostGIS uses GiST (Generalized Search Tree), which stores a bounding box per row in a tree, so "which boxes overlap this 10 km square" is answered by walking a few tree pages.
• geography vs geometry: geometry treats coordinates as flat x/y, so distances come out in degrees. geography measures on the Earth's surface, so ST_DWithin(a, b, 5000) means 5,000 metres. With raw lat, lon input and a kilometre requirement, geography is the correct type.
• SRID 4326 is the identifier for plain WGS 84 (World Geodetic System 1984) longitude/latitude, the coordinate system GPS reports.
• A lateral join (LATERAL) is a subquery that may reference columns of the rows to its left, so it runs once per search row. That is exactly the "top N per group" shape.
• KNN (k-nearest-neighbour) ordering: ORDER BY a <-> b LIMIT k lets PostGIS return rows nearest-first straight from the index.

Design decisions.

Runnable solution (PostgreSQL 16, PostGIS 3.4)

The script creates the tables, generates deterministic synthetic data with a pinned seed, builds the indexes, runs both queries and checks the plan. Run it with psql -f in a database where PostGIS is installed.

-- Setup: tables as given, plus a surrogate search_id and a city boundary table.
CREATE EXTENSION IF NOT EXISTS postgis;
DROP TABLE IF EXISTS listings, searches, cities;

CREATE TABLE listings (
  listing_id bigint PRIMARY KEY,
  lat double precision NOT NULL,
  lon double precision NOT NULL,
  is_active boolean NOT NULL,
  price_usd numeric(10,2) NOT NULL,
  geog geography(Point,4326)
    GENERATED ALWAYS AS (ST_SetSRID(ST_MakePoint(lon, lat), 4326)::geography) STORED
);
CREATE TABLE searches (
  search_id bigserial PRIMARY KEY,
  lat double precision NOT NULL,
  lon double precision NOT NULL,
  search_time timestamptz NOT NULL
);
CREATE TABLE cities (city text PRIMARY KEY, geom geometry(Polygon,4326) NOT NULL);

-- Deterministic synthetic data (PostgreSQL 16).
SELECT setseed(0.42);
INSERT INTO cities VALUES
 ('Austin', ST_MakeEnvelope(-97.90, 30.15, -97.60, 30.45, 4326)),
 ('Denver', ST_MakeEnvelope(-105.10, 39.60, -104.80, 39.90, 4326));
-- Austin: 3000 listings, Denver: 60 listings, ~80% active.
INSERT INTO listings (listing_id, lat, lon, is_active, price_usd)
SELECT g, 30.15 + random()*0.30, -97.90 + random()*0.30, random() < 0.8, round((50 + random()*450)::numeric, 2)
FROM generate_series(1, 3000) g;
INSERT INTO listings (listing_id, lat, lon, is_active, price_usd)
SELECT g, 39.60 + random()*0.30, -105.10 + random()*0.30, random() < 0.8, round((50 + random()*450)::numeric, 2)
FROM generate_series(3001, 3060) g;
-- 500 searches per city.
INSERT INTO searches (lat, lon, search_time)
SELECT 30.15 + random()*0.30, -97.90 + random()*0.30, timestamptz '2026-01-01' + g * interval '1 minute'
FROM generate_series(1, 500) g;
INSERT INTO searches (lat, lon, search_time)
SELECT 39.60 + random()*0.30, -105.10 + random()*0.30, timestamptz '2026-01-01' + g * interval '1 minute'
FROM generate_series(1, 500) g;

-- Indexes: partial GiST on active listings, GiST on city polygons.
CREATE INDEX listings_active_geog_gix ON listings USING gist (geog) WHERE is_active;
CREATE INDEX cities_geom_gix ON cities USING gist (geom);
ANALYZE listings; ANALYZE searches; ANALYZE cities;

-- Query 1: up to 5 nearest active listings within 5 km, per search.
CREATE TEMP TABLE search_results AS
SELECT s.search_id, s.search_time, COALESCE(c.city, 'unassigned') AS city,
       nn.listing_id, nn.distance_km, nn.price_usd
FROM searches s
CROSS JOIN LATERAL (SELECT ST_SetSRID(ST_MakePoint(s.lon, s.lat), 4326) AS pt) p
LEFT JOIN cities c ON ST_Covers(c.geom, p.pt)
LEFT JOIN LATERAL (
  SELECT l.listing_id, l.price_usd,
         round((ST_Distance(l.geog, p.pt::geography) / 1000.0)::numeric, 3) AS distance_km
  FROM listings l
  WHERE l.is_active
    AND ST_DWithin(l.geog, p.pt::geography, 5000)
  ORDER BY l.geog <-> p.pt::geography
  LIMIT 5
) nn ON true;

SELECT search_id, city, listing_id, distance_km, price_usd
FROM search_results WHERE search_id IN (1, 501)
ORDER BY search_id, distance_km;

-- Query 2: % of searches per city with at least 3 results.
SELECT city,
       count(*) AS searches,
       count(*) FILTER (WHERE n_results >= 3) AS with_3plus,
       round(100.0 * count(*) FILTER (WHERE n_results >= 3) / count(*), 1) AS pct_3plus
FROM (SELECT search_id, city, count(listing_id) AS n_results
      FROM search_results GROUP BY search_id, city) per_search
GROUP BY city ORDER BY city;

-- Plan check for the lateral probe.
EXPLAIN (COSTS OFF)
SELECT l.listing_id FROM listings l
WHERE l.is_active AND ST_DWithin(l.geog, ST_SetSRID(ST_MakePoint(-97.75, 30.30), 4326)::geography, 5000)
ORDER BY l.geog <-> ST_SetSRID(ST_MakePoint(-97.75, 30.30), 4326)::geography
LIMIT 5;

Output of that script (the empty setseed result and the NOTICE lines from the DROP ... IF EXISTS are omitted):

 search_id |  city  | listing_id | distance_km | price_usd
-----------+--------+------------+-------------+-----------
         1 | Austin |       1985 |       0.117 |    311.20
         1 | Austin |        925 |       0.197 |    239.97
         1 | Austin |       1743 |       0.398 |    223.44
         1 | Austin |       2100 |       0.554 |    445.05
         1 | Austin |        267 |       0.627 |    129.20
       501 | Denver |       3034 |       0.776 |    347.95
       501 | Denver |       3041 |       1.057 |    278.69
       501 | Denver |       3010 |       2.690 |    306.64
       501 | Denver |       3036 |       3.902 |    432.48

  city  | searches | with_3plus | pct_3plus
--------+----------+------------+-----------
 Austin |      500 |        500 |     100.0
 Denver |      500 |        373 |      74.6

 Limit
   ->  Sort
         Sort Key: ((geog <-> '...'::geography))
         ->  Bitmap Heap Scan on listings l
               Recheck Cond: is_active
               Filter: st_dwithin(geog, '...'::geography, '5000'::double precision, true)
               ->  Bitmap Index Scan on listings_active_geog_gix
                     Index Cond: (geog && _st_expand('...'::geography, '5000'::double precision))

(The long hex point literals in the plan are abbreviated to '...'.)

Reading the result

• Search 501 in Denver returned only 4 rows, which is correct: only 4 active listings sit within 5 km. The LIMIT 5 is a ceiling, not a guarantee.
• The density difference drives the percentage. Austin has 3,000 listings in a 0.30 by 0.30 degree box, so every Austin search clears 3 results. Denver has only 60 listings in a box of the same size, 46 of them active in this run: spread over that box, which at Denver's latitude works out to about 25.6 km by 33.4 km, roughly 856 km², that is about 0.054 active listings per km². A 5 km search radius is a circle of about pi x 5^2 = 78.5 km², so a search is expected to find about 0.054 x 78.5, roughly 4.2 active listings, close to the 3-result bar the query checks against. With the expected count that close to the bar, whether any one search actually clears it varies a lot from search to search, and 373 of 500 Denver searches do: 373 / 500 = 74.6%.
• Reading the plan, line by line. Bitmap Index Scan on listings_active_geog_gix walks the partial GiST index and builds a bitmap (an in-memory list of candidate row locations) instead of visiting table rows directly; Index Cond: (geog && _st_expand(...)) is the actual index lookup, where _st_expand grows the search point into a bounding box covering the 5,000 m radius and && (PostGIS's "bounding boxes overlap" operator) is the cheap, index-supported test the GiST tree can answer. Bitmap Heap Scan on listings l then visits just those candidate rows in the table's own storage (the heap, as opposed to the index) to check the conditions the bitmap step could only approximate: Recheck Cond: is_active re-checks the real condition against each row once it gets there, because a bitmap only marks which pages might contain a match, and Filter: st_dwithin(...) re-checks the true distance the same way, since the bounding box only proves "close enough to be worth checking", not "actually within 5 km". The Bitmap Index Scan line is what proves the partial index was actually used; the Recheck/Filter lines are why the plan is still correct even though its first pass over the box is approximate.
• The plan confirms the partial index is used: Bitmap Index Scan on listings_active_geog_gix with the && bounding-box condition that ST_DWithin generated. The planner chose a bitmap scan plus a small sort here rather than a pure KNN index walk, because a 5 km radius holds few enough rows that sorting them is cheap. Either plan touches only nearby rows, which is the point.

Complexity

Let S be the number of searches, L the number of active listings and m the listings inside one 5 km radius. Each probe costs roughly O(log L + m log m) (tree descent, then sort the candidates), so the whole query is about O(S (log L + m log m)); with m capped by the LIMIT 5 in practice this stays cheap, but the honest bound keeps the log m term rather than dropping it. The naive cross join is O(S * L): with 1,000 searches and 3,000 listings that is 3,000,000 distance computations, and it grows with both tables.

Edge cases and pitfalls

• Longitude first. ST_MakePoint(x, y) takes (lon, lat). Swapping them puts Austin in Antarctica and the query silently returns nothing.
• Distance units. ST_Distance on geography returns metres (divide by 1,000 for km). On geometry in SRID 4326 it returns degrees, which is meaningless as a distance.
• Sphere vs spheroid. <-> on geography orders by sphere distance while ST_Distance defaults to the more exact spheroid. The two can disagree by a fraction of a percent, which only matters for near-ties in ordering.
• Ties. Add , l.listing_id to the ORDER BY if results must be stable across runs when two listings are equidistant.
• Searches with no city. Counting only INNER JOIN matches hides searches outside the polygon table; the COALESCE(..., 'unassigned') bucket keeps them visible.
• search_time is unused by the spec, but the same pattern takes a time filter trivially (for example restricting to the last 7 days before the lateral step), which is how you would run it incrementally over a large searches table.
• At larger scale, run the lateral step in batches keyed by search_id or by day, and consider physically reordering the table by location so nearby listings share heap pages (the fixed-size blocks of storage that hold a table's actual rows, the "heap"; reading fewer distinct pages means fewer disk or cache reads). CLUSTER (the command that physically rewrites a table's heap in index order) needs a non-partial index, so it would use a separate full GiST index on geog.

Direct answer

Encode the rider's position at a precision whose cells are about the size of the search radius, collect the rider's cell plus enough rings of neighbouring cells to reach R metres in every direction, take every driver whose geohash starts with one of those cell strings as a candidate, then keep only candidates whose exact Haversine distance is at most R. The cells are a cheap coarse filter that may include extra drivers; the distance check makes the answer exact. The ring count must be computed from the cell's size in metres at that latitude, because a geohash cell's width in metres shrinks as you move away from the equator.

Terms used

• Geohash: a string that names a rectangular cell of the Earth. Each character adds 5 bits, alternately splitting longitude and latitude in half, so longer strings mean smaller cells, and a cell's children all start with the parent's string.
• Neighbour cells: the 8 cells touching a cell. A "ring" is one layer of neighbours around the block you already have.
• Haversine: the formula for the great-circle distance (shortest path over the surface of a sphere) between two latitude/longitude points.
• Candidate: a driver that passed the cheap cell filter but has not yet passed the exact distance check.

Approach

1. Cell size at a precision. With 5p bits, longitude gets ceil(5p / 2) bits and latitude floor(5p / 2). Cell height = 180 / 2^(latitude bits) degrees and width = 360 / 2^(longitude bits) degrees. Height in metres is height x 111,195 (metres per degree of latitude on a 6,371 km sphere); width in metres also multiplies by cos(latitude).
2. Rings needed. k_lat = ceil(R / height_m) and k_lon = ceil(R / width_m). The rider is somewhere inside the centre cell, so k cells in each direction always reach at least R beyond it. Width is computed at the most pole-ward latitude the search can reach, where cells are narrowest.
3. Enumerate cells by stepping from the centre cell's centre by whole cell sizes and re-encoding. Longitude is wrapped into [-180, 180) so the search crosses the antimeridian correctly; steps past a pole are skipped.
4. Candidates = drivers whose geohash prefix is in that cell set (in production, one lookup per cell in a map from cell string to driver IDs).
5. Exact filter with Haversine, then sort by distance.

Code

import math

BASE32 = "0123456789bcdefghjkmnpqrstuvwxyz"
EARTH_RADIUS_M = 6_371_008.8


def encode(lat, lon, precision):
    """Standard geohash: interleave longitude and latitude bits (longitude first), 5 bits per char."""
    lat_lo, lat_hi, lon_lo, lon_hi = -90.0, 90.0, -180.0, 180.0
    out, bits, ch, even = [], 0, 0, True
    while len(out) < precision:
        if even:
            mid = (lon_lo + lon_hi) / 2
            ch = (ch << 1) | (lon >= mid)
            lon_lo, lon_hi = (mid, lon_hi) if lon >= mid else (lon_lo, mid)
        else:
            mid = (lat_lo + lat_hi) / 2
            ch = (ch << 1) | (lat >= mid)
            lat_lo, lat_hi = (mid, lat_hi) if lat >= mid else (lat_lo, mid)
        even = not even
        bits += 1
        if bits == 5:
            out.append(BASE32[ch])
            bits, ch = 0, 0
    return "".join(out)


def cell_size_deg(precision):
    """Height and width of a cell in degrees. 5 bits per char, longitude gets the extra bit when odd."""
    total = 5 * precision
    lon_bits = (total + 1) // 2
    lat_bits = total // 2
    return 180.0 / 2 ** lat_bits, 360.0 / 2 ** lon_bits


def decode_center(gh):
    lat_lo, lat_hi, lon_lo, lon_hi = -90.0, 90.0, -180.0, 180.0
    even = True
    for c in gh:
        v = BASE32.index(c)
        for shift in range(4, -1, -1):
            bit = (v >> shift) & 1
            if even:
                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)
            even = not even
    return (lat_lo + lat_hi) / 2, (lon_lo + lon_hi) / 2


def haversine_m(lat1, lon1, lat2, lon2):
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dp, dl = p2 - p1, math.radians(lon2 - lon1)
    a = math.sin(dp / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dl / 2) ** 2
    return 2 * EARTH_RADIUS_M * math.asin(math.sqrt(min(1.0, a)))


def cells_covering(lat, lon, radius_m, precision):
    """Rider's cell plus k rings of neighbours, k chosen so the block reaches radius_m in every direction."""
    h_deg, w_deg = cell_size_deg(precision)
    h_m = h_deg * 111_195.0                       # metres per degree of latitude (pi * R / 180)
    # a cell's narrowest width is at its pole-ward edge, so size the rings for the worst latitude in reach
    worst_lat = min(89.9, abs(lat) + radius_m / 111_195.0 + h_deg)
    w_m = w_deg * 111_195.0 * math.cos(math.radians(worst_lat))
    k_lat = math.ceil(radius_m / h_m)
    k_lon = math.ceil(radius_m / w_m)
    c_lat, c_lon = decode_center(encode(lat, lon, precision))
    cells = set()
    for dy in range(-k_lat, k_lat + 1):
        nlat = c_lat + dy * h_deg
        if nlat <= -90 or nlat >= 90:             # past a pole: no cell there
            continue
        for dx in range(-k_lon, k_lon + 1):
            nlon = (c_lon + dx * w_deg + 180.0) % 360.0 - 180.0   # wrap across the antimeridian
            cells.add(encode(nlat, nlon, precision))
    return cells, (k_lat, k_lon), (h_m, w_m)


def nearby_drivers(drivers, rider_lat, rider_lon, radius_m, precision=7):
    cells, rings, size = cells_covering(rider_lat, rider_lon, radius_m, precision)
    # in production this is one lookup per cell in a map of cell -> driver ids; a scan stands in for it here
    candidates = [d for d in drivers if d["geohash"][:precision] in cells]
    hits = []
    for d in candidates:
        dist = haversine_m(rider_lat, rider_lon, d["lat"], d["lon"])
        if dist <= radius_m:
            hits.append((round(dist, 1), d["id"]))
    return sorted(hits), len(cells), rings, size, len(candidates)


# A few blocks around Market St and 4th St, San Francisco
rider = (37.78486, -122.40540)
raw = [
    ("d1", 37.78520, -122.40480),   # same block
    ("d2", 37.78610, -122.40710),   # across Market
    ("d3", 37.78300, -122.40300),   # two blocks south-east
    ("d4", 37.78790, -122.40140),   # Union Square-ish, just inside or outside?
    ("d5", 37.78050, -122.41000),   # further south-west
    ("d6", 37.79500, -122.39400),   # Embarcadero direction, well outside
]
drivers = [{"id": i, "lat": la, "lon": lo, "geohash": encode(la, lo, 9)} for i, la, lo in raw]

R = 400
print("rider geohash (p7):", encode(*rider, 7))
for prec in (6, 7):
    hits, n_cells, rings, (h_m, w_m), n_cand = nearby_drivers(drivers, *rider, R, precision=prec)
    print(f"p{prec}: cell ~{h_m:.0f} m tall x {w_m:.0f} m wide, rings (lat, lon) = {rings}, "
          f"cells = {n_cells}, candidates = {n_cand}, within {R} m = {hits}")
print("all distances:", {d['id']: round(haversine_m(*rider, d['lat'], d['lon']), 1) for d in drivers})

# Correctness check: brute force over 20,000 seeded random points around the rider
import random
rng = random.Random(7)
pts = [{"id": n, "lat": rider[0] + rng.uniform(-0.01, 0.01), "lon": rider[1] + rng.uniform(-0.012, 0.012)}
       for n in range(20_000)]
for p in pts:
    p["geohash"] = encode(p["lat"], p["lon"], 9)
got = {i for _, i in nearby_drivers(pts, *rider, R, precision=7)[0]}
truth = {p["id"] for p in pts if haversine_m(*rider, p["lat"], p["lon"]) <= R}
print("brute-force check:", len(got), "found,", len(truth), "true, missed =", len(truth - got))

# Antimeridian: a rider just west of 180 degrees must see a driver just east of it
fiji_rider = (-16.50, 179.9990)
fiji_driver = {"id": "east", "lat": -16.50, "lon": -179.9995, "geohash": encode(-16.50, -179.9995, 9)}
print("across antimeridian:", nearby_drivers([fiji_driver], *fiji_rider, 300, precision=7)[0],
      "| naive same-prefix?", encode(*fiji_rider, 7)[:3], "vs", fiji_driver["geohash"][:3])

Output:

rider geohash (p7): 9q8yywe
p6: cell ~611 m tall x 965 m wide, rings (lat, lon) = (1, 1), cells = 9, candidates = 5, within 400 m = [(64.9, 'd1'), (203.3, 'd2'), (295.4, 'd3')]
p7: cell ~153 m tall x 121 m wide, rings (lat, lon) = (3, 4), cells = 63, candidates = 4, within 400 m = [(64.9, 'd1'), (203.3, 'd2'), (295.4, 'd3')]
all distances: {'d1': 64.9, 'd2': 203.3, 'd3': 295.4, 'd4': 487.7, 'd5': 631.2, 'd6': 1508.2}
brute-force check: 2124 found, 2124 true, missed = 0
across antimeridian: [(159.9, 'east')] | naive same-prefix? rvp vs 2j0

Key points

• Both precisions return the same three drivers, as they must: the distance check makes the result exact, and precision only changes how much work the cell filter does.
• Precision 6 (about 611 m x 965 m cells) needs just 1 ring, 9 cells. The 9 cells cover 3 x 611 m by 3 x 965 m, about 5.3 km², versus a search circle of pi x 0.4² = 0.50 km², so they over-fetch roughly 10x the area.
• Precision 7 (about 153 m x 121 m) needs 3 rings north-south and 4 east-west, 63 cells, covering 7 x 153 m by 9 x 121 m, about 1.17 km², only about 2.3x the circle. Fewer wasted candidates, 7x more cell lookups. In a dense downtown with thousands of drivers per km², the finer precision wins; in a sparse suburb, the coarse one does.
• d4 at 487.7 m was a candidate at both precisions but was removed by the Haversine check, which is exactly the false positive the second stage exists for.
• The brute-force check over 20,000 seeded points found all 2,124 true points and missed none, which tests the ring arithmetic.
• Across the antimeridian, the rider's cell starts with rvp and the driver's with 2j0: no shared prefix at all, yet the driver is 159.9 m away. The longitude wrap in step 3 is what finds them.

Complexity

• Encoding: O(p) for precision p (7 characters = 35 bit steps).
• Cells enumerated: (2 k_lat + 1)(2 k_lon + 1), each O(p) to encode.
• With a cell-to-drivers hash map, candidate lookup is O(number of cells + candidates). The linear scan in the demo is O(n) and stands in for that map.
• Exact filtering: O(candidates), each a constant-cost Haversine.

Edge cases

• Rider near a cell edge: handled by the rings; never search the rider's cell alone.
• Antimeridian: handled by wrapping longitude; code that only compares prefixes misses these drivers.
• Poles: cells north of +90 or south of -90 do not exist and are skipped; the cosine term makes width tiny, so k_lon grows. Near the poles a sphere-native index (S2 or H3) is the better tool.
• Stale positions: a driver's stored geohash may lag their real position by one ping interval (for example 4 s at 15 m/s = 60 m); either add that to R for the cell filter or re-check against the freshest position.
• Straight-line distance is not travel distance: a driver 200 m away across a river may be 3 km by road. Real matching re-ranks candidates by estimated travel time.

Direct answer

Convert both points to radians, apply the haversine formula to get the central angle between them, and multiply by the Earth's mean radius (6,371,008.8 m). Validate that each input is a finite real number with latitude in [-90, 90] and longitude in [-180, 180]. The antimeridian (the 180° line) needs no special branch, because the formula uses the sine of half the longitude difference squared, which gives the same value for a 359.8° difference as for a 0.2° one. For the example, San Francisco to Oakland comes out at about 13,430 m.

The formula

Haversine distance is the great-circle distance: the shortest path between two points along the surface of a sphere. With latitude φ, longitude λ, both in radians, and radius R:

Here a is a number between 0 (same point) and 1 (opposite sides of the Earth), and 2·arcsin(√a) is the angle between the two points as seen from the Earth's centre.

Code

import math
from numbers import Real

EARTH_RADIUS_M = 6_371_008.8  # IUGG mean Earth radius, in meters (IUGG: the International Union
                               # of Geodesy and Geophysics, which publishes the reference figures
                               # for the Earth's size and shape)


def _check(name, value, low, high):
    # bool is a subclass of int in Python, so reject it explicitly
    if isinstance(value, bool) or not isinstance(value, Real):
        raise TypeError(f"{name} must be a real number, got {type(value).__name__}")
    value = float(value)
    if not math.isfinite(value):
        raise ValueError(f"{name} must be finite, got {value}")
    if not low <= value <= high:
        raise ValueError(f"{name}={value} is outside [{low}, {high}]")
    return value


def haversine_m(lat1, lon1, lat2, lon2, radius_m=EARTH_RADIUS_M):
    """Great-circle distance in meters between two (lat, lon) points in decimal degrees."""
    lat1 = _check("lat1", lat1, -90.0, 90.0)
    lat2 = _check("lat2", lat2, -90.0, 90.0)
    lon1 = _check("lon1", lon1, -180.0, 180.0)
    lon2 = _check("lon2", lon2, -180.0, 180.0)

    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    dphi = phi2 - phi1
    dlmb = math.radians(lon2 - lon1)  # may be ~360 deg across the antimeridian; sin^2 below is periodic

    a = math.sin(dphi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(dlmb / 2) ** 2
    a = min(1.0, max(0.0, a))  # float rounding can push a a hair outside [0, 1]
    return 2 * radius_m * math.asin(math.sqrt(a))


if __name__ == "__main__":
    cases = {
        "SF -> Oakland": (37.7749, -122.4194, 37.8044, -122.2712),
        "SF -> LA": (37.7749, -122.4194, 34.0522, -118.2437),
        "antimeridian (0,179.9) -> (0,-179.9)": (0.0, 179.9, 0.0, -179.9),
        "Fiji -> Samoa (crosses 180)": (-17.7134, 178.0650, -13.7590, -172.1046),
        "same point": (51.5, -0.12, 51.5, -0.12),
        "antipodal (0,0) -> (0,180)": (0.0, 0.0, 0.0, 180.0),
        "pole to pole": (90.0, 0.0, -90.0, 0.0),
    }
    for name, args in cases.items():
        d = haversine_m(*args)
        print(f"{name:40s} {d:14.1f} m  ({d / 1000:.3f} km)")

    for bad in [(91, 0, 0, 0), (0, 200, 0, 0), (float("nan"), 0, 0, 0), ("37.7", 0, 0, 0), (True, 0, 0, 0)]:
        try:
            haversine_m(*bad)
        except (TypeError, ValueError) as e:
            print(f"rejected {bad!r}: {type(e).__name__}: {e}")

Output:

SF -> Oakland                                   13429.6 m  (13.430 km)
SF -> LA                                       559121.3 m  (559.121 km)
antimeridian (0,179.9) -> (0,-179.9)            22239.0 m  (22.239 km)
Fiji -> Samoa (crosses 180)                   1139984.3 m  (1139.984 km)
same point                                          0.0 m  (0.000 km)
antipodal (0,0) -> (0,180)                   20015114.4 m  (20015.114 km)
pole to pole                                 20015114.4 m  (20015.114 km)
rejected (91, 0, 0, 0): ValueError: lat1=91.0 is outside [-90.0, 90.0]
rejected (0, 200, 0, 0): ValueError: lon1=200.0 is outside [-180.0, 180.0]
rejected (nan, 0, 0, 0): ValueError: lat1 must be finite, got nan
rejected ('37.7', 0, 0, 0): TypeError: lat1 must be a real number, got str
rejected (True, 0, 0, 0): TypeError: lat1 must be a real number, got bool

Key points

• Validation: bool is a subclass of int in Python, so True would otherwise pass as 1.0; strings are rejected rather than silently parsed; NaN and infinity are rejected because they would propagate into a NaN distance.
• Antimeridian: (0, 179.9) to (0, -179.9) has a raw longitude difference of -359.8°. Half of that is -179.9°, and sin²(-179.9°) = sin²(0.1°), exactly as for a 0.2° gap. The result is 22,239.0 m, which matches 0.2° of the equator: 0.2 × π / 180 × 6,371,008.8 ≈ 22,239 m. The Fiji-to-Samoa case also crosses 180° and returns about 1,140 km, not tens of thousands.
• Clamping a to [0, 1]: for nearly antipodal points floating-point rounding can make a a hair above 1, and asin of a value above 1 raises ValueError. The clamp makes the function total (defined for every valid input, with none of them causing it to error out or return an undefined result).
• Units: return metres by default. The kilometre version is the same function divided by 1,000; San Francisco (37.7749, -122.4194) to Los Angeles (34.0522, -118.2437) gives 559,121.3 m, or 559.1 km, a useful sanity check because the commonly quoted figure is about 559 km.
• Radius choice: 6,371,008.8 m is the mean Earth radius. Some code uses 6,371,000 m or 6,378,137 m (the equatorial radius). The last one inflates every distance by about 0.11%, so be consistent across services that compare distances.

Complexity

O(1) time and O(1) memory per call: a fixed number of trig operations. For a million pairs, vectorise the same formula with NumPy arrays rather than looping in Python.

Edge cases

How accurate is it? (Sphere vs the real Earth)

The Earth is slightly flattened (an ellipsoid). The standard ellipsoid model used by GPS is WGS84. Comparing against an exact ellipsoidal geodesic distance (a geodesic is the shortest path between two points across a curved surface, here the ellipsoid rather than a sphere) using the pyproj library (pip install pyproj, and it imports the function above from haversine.py):

import math
from pyproj import Geod
from haversine import haversine_m

geod = Geod(ellps="WGS84")
for name, (la1, lo1, la2, lo2) in {
    "SF -> Oakland": (37.7749, -122.4194, 37.8044, -122.2712),
    "SF -> LA": (37.7749, -122.4194, 34.0522, -118.2437),
    "antimeridian": (0.0, 179.9, 0.0, -179.9),
    "pole to pole": (90.0, 0.0, -90.0, 0.0),
}.items():
    _, _, g = geod.inv(lo1, la1, lo2, la2)  # pyproj takes lon, lat order
    h = haversine_m(la1, lo1, la2, lo2)
    print(f"{name:14s} haversine={h:12.1f} m  WGS84={g:12.1f} m  diff={100 * (h - g) / g:+.3f}%")

# Pitfall: equirectangular approximation without wrapping delta-longitude
def equirect_naive_m(lat1, lon1, lat2, lon2, R=6_371_008.8):
    x = math.radians(lon2 - lon1) * math.cos(math.radians((lat1 + lat2) / 2))
    y = math.radians(lat2 - lat1)
    return R * math.hypot(x, y)
print(f"naive equirectangular across 180: {equirect_naive_m(0.0, 179.9, 0.0, -179.9) / 1000:.1f} km")

Output:

SF -> Oakland  haversine=     13429.6 m  WGS84=     13458.2 m  diff=-0.212%
SF -> LA       haversine=    559121.3 m  WGS84=    559042.3 m  diff=+0.014%
antimeridian   haversine=     22239.0 m  WGS84=     22263.9 m  diff=-0.112%
pole to pole   haversine=  20015114.4 m  WGS84=  20003931.5 m  diff=+0.056%
naive equirectangular across 180: 40008.0 km

The error is a few tenths of a percent at most: 28.6 m on the 13.4 km San Francisco to Oakland trip. That is fine for "which drivers are nearby" filtering and for display. It is not fine for surveying, legal boundaries or billing by distance, where you would use an ellipsoidal geodesic (Karney's algorithm, a modern method that always converges, including for antipodal, nearly opposite, points, and which pyproj.Geod implements, or Vincenty's, which can fail to converge for those same points).

Trade-offs and pitfalls

• Naive flat formulas break at 180°. The equirectangular shortcut above, without wrapping the longitude difference, reports 40,008 km for two points 22 km apart. If you use a flat approximation for speed inside a city, wrap Δλ into [-180°, 180°] first.
• Swapped arguments are the most common real bug. Many libraries (GeoJSON, PostGIS, pyproj) order coordinates as (lon, lat). Keyword arguments or a named tuple protect you.
• Straight-line distance is not travel distance. For ride-hailing ETAs, haversine only shortlists candidates; the road network decides.