Geospatial and Time-Series Data Questions

Direct answer

For 100k drivers in one process, build a scikit-learn BallTree over driver positions using the haversine metric (great-circle distance on a sphere), converted to radians with latitude first, and call query(k=...). Build costs O(n log n), a query is roughly O(log n + k) on spread-out 2-D data, and memory is O(n), a few MB. Because these trees are static, handle moving drivers by keeping live positions in a plain array, rebuilding the tree as a snapshot every couple of seconds, and over-fetching candidates from the snapshot so they can be re-ranked by where each driver is NOW.

Terms used

• k-nearest-neighbour (kNN) query: "give me the k points closest to this one."
• Haversine distance: the shortest distance between two points along the surface of a sphere (the great-circle distance). It needs angles in radians.
• KDTree: a tree that repeatedly splits space with axis-aligned cuts (first by x, then by y, and so on). It assumes a flat, Euclidean space.
• BallTree: a tree that groups points into nested "balls" (a centre and a radius). It only needs a valid distance function, so it works with haversine.
• Projection: converting (lat, lon) degrees into flat x/y metres so ordinary Euclidean distance is meaningful.

Approach

1. Pick coordinates that make the distance honest. Degrees of longitude shrink with latitude: at 37.75° N one degree of longitude is about 111.2 × cos(37.75°) ≈ 87.9 km, while one degree of latitude is about 111.2 km. Feeding raw (lat, lon) degrees into a Euclidean KDTree therefore ranks neighbours with a stretched ruler. Three correct options exist; this answer commits to the first (BallTree + haversine) and uses it for the rest of the code below, the other two are useful context, not the default:
   • BallTree + haversine on [lat, lon] in radians. Exact on the sphere, works everywhere, including across the 180° meridian (the antimeridian).
   • KDTree on projected metres (for example a local equirectangular projection centred on the city: x = R·lon·cos(lat0), y = R·lat). Fast and fine for one metro area, but distortion grows as you move away from the centre latitude, and it breaks at the antimeridian.
   • KDTree on 3-D unit vectors (x, y, z on the unit sphere), a niche trick for when only a plain Euclidean KDTree library is available. Straight-line chord distance (a straight segment connecting two points, cutting through the Earth's interior rather than following its curved surface) through the Earth grows monotonically (it only ever increases, never decreases, as the true great-circle distance increases) with great-circle distance, so the neighbour ORDER is exact globally. Convert back with angle = 2·asin(chord/2).
2. Scale the result. BallTree(metric="haversine") returns distances in radians on a unit sphere; multiply by the Earth radius (6,371,008.8 m, the mean radius) to get metres.
3. Verify against brute force. For a new index, always compare with an O(n) exact scan on a sample. It catches the classic bug of passing [lon, lat] or degrees.

Code: static kNN and the coordinate choices

Requires numpy and scikit-learn.

import numpy as np
from sklearn.neighbors import BallTree, KDTree

EARTH_RADIUS_M = 6_371_008.8
rng = np.random.default_rng(42)

# 100k synthetic drivers in a San Francisco Bay Area box, columns = [lat, lon] in degrees
N = 100_000
drivers_deg = np.column_stack([
    rng.uniform(37.60, 37.90, N),
    rng.uniform(-122.55, -122.20, N),
])


def build_index(latlon_deg):
    # haversine metric in scikit-learn expects [lat, lon] in RADIANS, latitude first
    # leaf_size caps how many points sit in a leaf before the tree stops splitting
    # and falls back to a linear scan inside it; it trades memory for query speed
    # (memory drops roughly in proportion to leaf_size) without changing which
    # neighbours are returned. 40 is scikit-learn's own default.
    return BallTree(np.radians(latlon_deg), metric="haversine", leaf_size=40)


def k_nearest(tree, rider_lat, rider_lon, k=5):
    q = np.radians([[rider_lat, rider_lon]])
    dist_rad, idx = tree.query(q, k=k)          # distances come back in radians (unit sphere)
    return idx[0], dist_rad[0] * EARTH_RADIUS_M  # scale by Earth radius to get meters


def brute_force(latlon_deg, rider_lat, rider_lon, k=5):
    lat, lon = np.radians(latlon_deg[:, 0]), np.radians(latlon_deg[:, 1])
    p, l = np.radians(rider_lat), np.radians(rider_lon)
    a = np.sin((lat - p) / 2) ** 2 + np.cos(p) * np.cos(lat) * np.sin((lon - l) / 2) ** 2
    d = 2 * EARTH_RADIUS_M * np.arcsin(np.sqrt(np.clip(a, 0, 1)))
    order = np.argsort(d)[:k]
    return order, d[order]


tree = build_index(drivers_deg)
rider = (37.7749, -122.4194)
idx, dist_m = k_nearest(tree, *rider, k=5)
bf_idx, bf_d = brute_force(drivers_deg, *rider, k=5)
print("BallTree   ids:", idx.tolist())
print("BallTree   m  :", np.round(dist_m, 1).tolist())
print("brute force ids match:", idx.tolist() == bf_idx.tolist(), "| max |diff| m:", float(np.max(np.abs(dist_m - bf_d))))

# KDTree cannot use haversine; it only supports Minkowski-family metrics
# (distance formulas built from raw coordinate differences, such as Euclidean
# or Manhattan distance, that assume a flat space)
print("'haversine' in KDTree.valid_metrics:", "haversine" in KDTree.valid_metrics)

# Alternative 1: local equirectangular projection to meters, then KDTree (Euclidean)
lat0 = np.radians(rider[0])
def project(latlon_deg):
    lat, lon = np.radians(latlon_deg[:, 0]), np.radians(latlon_deg[:, 1])
    return np.column_stack([EARTH_RADIUS_M * lon * np.cos(lat0), EARTH_RADIUS_M * lat])
kd = KDTree(project(drivers_deg))
kd_d, kd_i = kd.query(project(np.array([rider])), k=5)
print("KDTree(projected) ids match:", kd_i[0].tolist() == bf_idx.tolist(), "| max |diff| m:", round(float(np.max(np.abs(kd_d[0] - bf_d))), 3))

# Alternative 2: 3D unit vectors + KDTree. Chord length is monotonic in great-circle angle,
# so the neighbor ORDER is exact anywhere on Earth (poles, antimeridian included).
def to_xyz(latlon_deg):
    lat, lon = np.radians(latlon_deg[:, 0]), np.radians(latlon_deg[:, 1])
    return np.column_stack([np.cos(lat) * np.cos(lon), np.cos(lat) * np.sin(lon), np.sin(lat)])
kd3 = KDTree(to_xyz(drivers_deg))
chord, i3 = kd3.query(to_xyz(np.array([rider])), k=5)
arc_m = 2 * np.arcsin(chord[0] / 2) * EARTH_RADIUS_M
print("KDTree(xyz) ids match:", i3[0].tolist() == bf_idx.tolist(), "| max |diff| m:", round(float(np.max(np.abs(arc_m - bf_d))), 6))

# Antimeridian check: two drivers either side of 180 degrees, rider at lon 179.99
edge = np.radians([[0.0, 179.995], [0.0, -179.995], [0.0, 179.0]])
edge_tree = BallTree(edge, metric="haversine")
d, i = edge_tree.query(np.radians([[0.0, 179.99]]), k=3)
print("antimeridian order:", i[0].tolist(), "meters:", np.round(d[0] * EARTH_RADIUS_M, 1).tolist())

Output:

BallTree   ids: [49604, 87224, 30155, 78708, 38563]
BallTree   m  : [37.4, 48.0, 72.1, 94.9, 106.4]
brute force ids match: True | max |diff| m: 0.0
'haversine' in KDTree.valid_metrics: False
KDTree(projected) ids match: True | max |diff| m: 0.0
KDTree(xyz) ids match: True | max |diff| m: 0.0
antimeridian order: [0, 1, 2] meters: [556.0, 1667.9, 110083.1]

All three correct methods agree with the brute-force scan, and KDTree confirms it has no haversine metric. The antimeridian test puts the driver at longitude -179.995 (1,667.9 m away) ahead of the one at 179.0 (110 km away), which a naive "subtract the longitudes" approach gets backwards.

Supporting dynamic updates

scikit-learn's BallTree/KDTree and SciPy's KDTree have no insert, delete or move operation: changing a point means rebuilding. With drivers sending a GPS ping every 4 s, 100,000 drivers produce 100,000 / 4 = 25,000 updates per second, so updating the tree per ping is not an option. The design I would commit to for a single-city, in-process service:

• Live state is an array, one row per driver. An update is an O(1) write; availability (online, not on a trip) is a boolean per driver.
• The tree is a snapshot rebuilt on a timer (every 2 s here) off the request path, then swapped in atomically, so a reader never sees a half-built tree.
• Queries over-fetch and re-rank. Ask the snapshot for k + 20 candidates, drop unavailable drivers, then rank those candidates by haversine distance from their CURRENT positions.

import threading
import numpy as np
from sklearn.neighbors import BallTree

EARTH_RADIUS_M = 6_371_008.8


def haversine_m(lat1, lon1, lat2, lon2):
    lat1, lon1, lat2, lon2 = map(np.radians, (lat1, lon1, lat2, lon2))
    a = np.sin((lat2 - lat1) / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin((lon2 - lon1) / 2) ** 2
    return 2 * EARTH_RADIUS_M * np.arcsin(np.sqrt(np.clip(a, 0, 1)))


class DriverIndex:
    """Live positions are updated in O(1); the BallTree is a periodic snapshot.

    Queries over-fetch from the snapshot, then re-rank the candidates by their
    CURRENT position, so a moved driver is scored where it is now.
    """

    def __init__(self, latlon_deg):
        self.pos = latlon_deg.copy()                 # live [lat, lon] per driver id
        self.available = np.ones(len(latlon_deg), bool)
        self._lock = threading.Lock()
        self.rebuild()

    def update(self, driver_id, lat, lon):           # O(1): no tree touched
        self.pos[driver_id] = (lat, lon)

    def set_available(self, driver_id, flag):
        self.available[driver_id] = flag

    def rebuild(self):                               # O(n log n), run on a timer off the hot path
                                                       # (the latency-sensitive code that handles each
                                                       # incoming request; kept separate so a query is
                                                       # never blocked on this rebuild)
        snapshot = self.pos.copy()
        tree = BallTree(np.radians(snapshot), metric="haversine")
        with self._lock:                             # atomic swap (the tree reference changes in one
                                              # indivisible step): readers never see a half-built tree
            self._tree = tree

    def k_nearest(self, lat, lon, k=5, overfetch=20):
        with self._lock:
            tree = self._tree
        _, cand = tree.query(np.radians([[lat, lon]]), k=k + overfetch)
        cand = cand[0][self.available[cand[0]]]      # drop drivers who went offline or got a trip
        live = self.pos[cand]
        d = haversine_m(lat, lon, live[:, 0], live[:, 1])
        order = np.argsort(d)[:k]
        return cand[order], d[order]


rng = np.random.default_rng(7)
N = 100_000
drivers = np.column_stack([rng.uniform(37.60, 37.90, N), rng.uniform(-122.55, -122.20, N)])
index = DriverIndex(drivers)

# Between two rebuilds every driver moves up to 30 m (15 m/s for 2 s) in a random direction
step_m = rng.uniform(0, 30, N)
heading = rng.uniform(0, 2 * np.pi, N)
dlat = np.degrees(step_m * np.cos(heading) / EARTH_RADIUS_M)
dlon = np.degrees(step_m * np.sin(heading) / (EARTH_RADIUS_M * np.cos(np.radians(drivers[:, 0]))))
for i in range(N):
    index.update(i, drivers[i, 0] + dlat[i], drivers[i, 1] + dlon[i])
index.set_available(123, False)

riders = np.column_stack([rng.uniform(37.65, 37.85, 1000), rng.uniform(-122.50, -122.25, 1000)])
for overfetch in (0, 20):
    exact = 0
    for lat, lon in riders:
        got, _ = index.k_nearest(lat, lon, k=5, overfetch=overfetch)
        d = haversine_m(lat, lon, index.pos[:, 0], index.pos[:, 1])
        d[~index.available] = np.inf
        truth = np.argsort(d)[:5]
        exact += got.tolist() == truth.tolist()
    print(f"stale snapshot, overfetch={overfetch:2d}: exact top-5 for {exact}/1000 riders")

index.rebuild()
exact = sum(
    index.k_nearest(lat, lon, 5)[0].tolist()
    == np.argsort(np.where(index.available, haversine_m(lat, lon, index.pos[:, 0], index.pos[:, 1]), np.inf))[:5].tolist()
    for lat, lon in riders
)
print(f"fresh snapshot, overfetch=20: exact top-5 for {exact}/1000 riders")

Output:

stale snapshot, overfetch= 0: exact top-5 for 586/1000 riders
stale snapshot, overfetch=20: exact top-5 for 1000/1000 riders
fresh snapshot, overfetch=20: exact top-5 for 1000/1000 riders

Why 20 extra candidates is enough here (and when it is not)

The simulated box is 0.30° × 111.2 km ≈ 33.4 km tall by 0.35° × 87.9 km ≈ 30.8 km wide, about 1,027 km², so density ρ ≈ 100,000 / 1,027 ≈ 97 drivers per km². For randomly scattered points, a circle of radius r centred on the query point holds on average rho x pi x r^2 drivers (density times the circle's area); setting that equal to k and solving for r gives the distance to the k-th nearest neighbour, roughly

which gives about 128 m for k = 5 and about 286 m for k = 25. A driver can only be wrongly left out if it sat beyond the 25th candidate at snapshot time and moved inside the true top 5 before the next rebuild. The gap is about 286 − 128 ≈ 158 m, and in 2 s at 15 m/s (54 km/h) a driver moves at most 30 m, so the over-fetch margin is roughly five times what motion can consume. Without over-fetch, only 586 of 1,000 riders got the exact top 5. If drivers are sparser (suburbs) or the rebuild interval is longer, re-derive the over-fetch from the same inequality rather than keeping 20.

Complexity

Edge cases

• Order and units: haversine in scikit-learn wants [lat, lon] in radians. Passing degrees or [lon, lat] returns wrong neighbours silently.
• Fewer than k drivers (or fewer than k + m after filtering): query raises if k exceeds n; clamp k to the number of points and return what exists.
• Antimeridian and poles: fine with haversine or 3-D vectors; broken with a flat projection.
• Ties and duplicates (drivers queued at an airport at the same coordinate): order among equals is arbitrary; add a tie-breaker such as time since last trip if fairness matters.
• Stale drivers: a driver whose last ping is older than, say, 3 missed intervals should be marked unavailable, not matched at an old position.

Trade-offs and pitfalls

• Straight-line distance is a candidate filter, not the answer. A driver 300 m away across a river or on the wrong side of a divided highway may be 10 minutes out. Production dispatch takes the kNN shortlist (say 20 to 50) and ranks it by road-network ETA (estimated time of arrival) from a routing engine.
• Projected KDTree vs BallTree: inside one city a projected KDTree gives the same neighbours (shown above) and is simpler to reason about in metres. Choose BallTree or 3-D vectors when the service spans many latitudes or the antimeridian (Fiji, or New Zealand together with its Chatham Islands).
• When to leave the in-process tree: once positions must be shared across many stateless API servers, or the fleet spans many cities, move the live index out of process. Common choices are a grid index (bucket drivers by geohash or H3 cell, where an update is an O(1) move between buckets and a query scans the rider's cell plus neighbouring cells), or Redis GEO commands (GEOADD to upsert a position, GEOSEARCH ... BYRADIUS ... ASC COUNT k to query). The snapshot-and-re-rank idea still applies to any index that is cheaper to rebuild than to update.
• Rebuild cost is O(n log n) per interval, so doubling the fleet or halving the interval roughly doubles background CPU. Measure it on your hardware before choosing the interval rather than assuming.

Direct answer

Planar (Euclidean) distance treats a small patch of Earth as flat; with the longitude difference scaled by cos(latitude) it is cheapest and, within a metropolitan radius (under about 50 km), agrees with Haversine to well under a metre. Haversine models Earth as a sphere; it is accurate to about 0.5% (the sphere-vs-ellipsoid difference), well-behaved everywhere, and a few trigonometric calls per pair. Vincenty models Earth as an ellipsoid (the WGS84 shape, slightly flattened at the poles) and is accurate to millimetres, but it iterates and can fail to converge for nearly opposite points (called antipodal points: think of two people standing on exactly opposite sides of the globe). For ride matching, use a planar approximation (or Haversine) to prefilter and rank, because road travel time, not geometric distance, decides the final choice, and keep Vincenty (or better, Karney's geodesic algorithm, a newer method that solves the same ellipsoid problem but always converges, including for antipodal points) for surveying-grade or billing-grade distances.

Terms used

• Great-circle distance: the shortest path between two points along the surface of a sphere.
• Ellipsoid: a squashed sphere; WGS84, the model GPS uses, has an equatorial radius of 6,378,137 m and a flattening of 1/298.257223563.
• Equirectangular (planar) approximation: x = delta-longitude x cos(mean latitude), y = delta-latitude, both in radians, then Pythagoras, times the Earth radius.
• Antimeridian: the 180 degree meridian where longitude jumps from +180 to -180.
• Vectorization: computing many distances in one array operation (NumPy, or SIMD, single instruction multiple data, CPU instructions) instead of a loop.

The three formulas

Haversine, with latitudes phi, longitudes lambda, and Earth radius R:

Planar approximation:

Vincenty solves the inverse geodesic problem (given two points, find the distance between them along the ellipsoid's surface, as opposed to the "direct" problem of finding an end point from a start, a bearing and a distance) on the ellipsoid by iterating on the longitude difference on an auxiliary sphere (a plain sphere Vincenty uses as scratch space, because the trigonometry is easier there; the result is then corrected back onto the ellipsoid) until it changes by less than a tolerance, usually within a few iterations for normal distances.

The code below is not meant to be reproduced in an interview; what matters is the shape of the method. U1 and U2 are the reduced latitude (each input latitude adjusted for the ellipsoid's flattening, so the auxiliary-sphere trigonometry comes out right). lam is the running estimate of the longitude difference on the auxiliary sphere, and the loop refines it until two successive estimates agree to within tol. Once lam has converged, u2, AA and BB are terms of a series expansion (an approximation built from a sum of shrinking correction terms) that turns the auxiliary-sphere angle into a distance on the true ellipsoid. The takeaway, not the algebra: iterate on longitude over a sphere, then correct with a series for the ellipsoid's actual shape.

Measured comparison

import math

R = 6_371_008.8                       # mean Earth radius, metres
A, F = 6_378_137.0, 1 / 298.257223563  # WGS84 ellipsoid
B = A * (1 - F)


def planar(lat1, lon1, lat2, lon2, wrap=True):
    """Equirectangular: treat a small patch as flat, shrink longitude by cos(mean latitude)."""
    dlon = lon2 - lon1
    if wrap:                                   # take the short way round across the antimeridian
        dlon = (dlon + 180.0) % 360.0 - 180.0
    x = math.radians(dlon) * math.cos(math.radians((lat1 + lat2) / 2))
    y = math.radians(lat2 - lat1)
    return R * math.hypot(x, y)


def haversine(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 * R * math.asin(math.sqrt(min(1.0, a)))


def vincenty(lat1, lon1, lat2, lon2, tol=1e-12, max_iter=200):
    """Vincenty inverse formula on the WGS84 ellipsoid. Returns (metres, iterations)."""
    L = math.radians(lon2 - lon1)
    U1 = math.atan((1 - F) * math.tan(math.radians(lat1)))
    U2 = math.atan((1 - F) * math.tan(math.radians(lat2)))
    sU1, cU1, sU2, cU2 = math.sin(U1), math.cos(U1), math.sin(U2), math.cos(U2)
    lam = L
    for i in range(1, max_iter + 1):
        sl, cl = math.sin(lam), math.cos(lam)
        s_sig = math.hypot(cU2 * sl, cU1 * sU2 - sU1 * cU2 * cl)
        if s_sig == 0:
            return 0.0, i
        c_sig = sU1 * sU2 + cU1 * cU2 * cl
        sig = math.atan2(s_sig, c_sig)
        s_alpha = cU1 * cU2 * sl / s_sig
        c2_alpha = 1 - s_alpha ** 2
        c2sm = c_sig - 2 * sU1 * sU2 / c2_alpha if c2_alpha else 0.0
        C = F / 16 * c2_alpha * (4 + F * (4 - 3 * c2_alpha))
        lam_prev = lam
        lam = L + (1 - C) * F * s_alpha * (sig + C * s_sig * (c2sm + C * c_sig * (-1 + 2 * c2sm ** 2)))
        if abs(lam - lam_prev) < tol:
            break
    else:
        return float("nan"), max_iter        # did not converge (near-antipodal points)
    u2 = c2_alpha * (A ** 2 - B ** 2) / B ** 2
    AA = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
    BB = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
    d_sig = BB * s_sig * (c2sm + BB / 4 * (c_sig * (-1 + 2 * c2sm ** 2)
                          - BB / 6 * c2sm * (-3 + 4 * s_sig ** 2) * (-3 + 4 * c2sm ** 2)))
    return B * AA * (sig - d_sig), i


cases = {
    "SF 1 km north-east":     (37.7749, -122.4194, 37.7813, -122.4113),
    "SF 10 km east":          (37.7749, -122.4194, 37.7749, -122.3058),
    "SF 50 km south-east":    (37.7749, -122.4194, 37.4600, -121.9500),
    "Oslo 50 km north":       (59.9139, 10.7522, 60.3630, 10.7522),
    "over the pole 4.5 km":       (89.98, 0.0, 89.98, 180.0),
    "antimeridian 2 km":      (-16.50, 179.99, -16.50, -179.9912),
}
print(f"{'case':22} {'vincenty m':>11} {'haversine':>10} {'planar':>9} {'planar-hav m':>12} {'iters':>5}")
for name, (a, b, c, d) in cases.items():
    v, it = vincenty(a, b, c, d)
    h, p = haversine(a, b, c, d), planar(a, b, c, d)
    print(f"{name:22} {v:11.1f} {100 * (h - v) / v:+9.3f}% {100 * (p - v) / v:+8.3f}% {p - h:+12.2f} {it:5d}")

# Ranking: does the cheap metric reorder candidates? 2,000 seeded drivers within ~8 km of SF
import random
rng = random.Random(3)
rider = (37.7749, -122.4194)
drivers = [(rider[0] + rng.uniform(-0.07, 0.07), rider[1] + rng.uniform(-0.09, 0.09)) for _ in range(2000)]
top = lambda fn: [i for _, i in sorted((fn(*rider, *d), i) for i, d in enumerate(drivers))[:10]]
print("top-10 identical, planar vs vincenty:", top(planar) == top(lambda *x: vincenty(*x)[0]))
print("top-10 identical, haversine vs vincenty:", top(haversine) == top(lambda *x: vincenty(*x)[0]))
# Naive degree-space Euclidean (no cos(lat) correction) for contrast
naive = lambda a, b, c, d: math.hypot(c - a, d - b)
print("top-10 identical, raw-degree euclidean vs vincenty:", top(naive) == top(lambda *x: vincenty(*x)[0]))
print("antimeridian, planar without the wrap:", round(planar(-16.50, 179.99, -16.50, -179.9912, wrap=False) / 1000, 1), "km")

Output:

case                    vincenty m  haversine    planar planar-hav m iters
SF 1 km north-east          1006.9    -0.028%   -0.028%        +0.00     4
SF 10 km east              10008.2    -0.237%   -0.237%        +0.00     4
SF 50 km south-east        54212.0    -0.061%   -0.061%        +0.14     4
Oslo 50 km north           50036.3    -0.197%   -0.197%        +0.00     1
over the pole 4.5 km        4467.8    -0.447%  +56.378%     +2538.79     1
antimeridian 2 km           2007.2    -0.139%   -0.139%        +0.00     4
top-10 identical, planar vs vincenty: True
top-10 identical, haversine vs vincenty: True
top-10 identical, raw-degree euclidean vs vincenty: False
antimeridian, planar without the wrap: 38379.8 km

What the numbers show

• Within a metro area, the formula choice barely matters for accuracy. Planar and Haversine differ by 0.00 to 0.14 m up to 54 km. Both differ from the ellipsoid (Vincenty) by 0.03% to 0.24%, a systematic error from modelling Earth as a sphere, at most 24 m on a 10 km trip. That error is nearly identical for every driver in the same area, so it does not change the ranking.
• Ranking: across 2,000 seeded drivers, the top 10 by planar and by Haversine are identical to the top 10 by Vincenty. Raw-degree Euclidean distance (no cos(latitude) factor) produces a different top 10, because at 37.8 degrees a degree of longitude is only about 79% as long as a degree of latitude. That is the planar mistake that actually matters.
• Poles: two points 4.5 km apart on opposite sides of the North Pole give +56% error with the planar formula, because a flat patch cannot represent paths that cross the pole. Haversine stays at -0.45%.
• Antimeridian: the planar formula without wrapping the longitude difference returns 38,380 km for a 2 km trip. Wrapping the longitude difference into [-180, 180) fixes it; Haversine is immune because sin²(delta-lambda / 2) is periodic.
• Numerical stability: Haversine is well-conditioned (small changes in the input only ever produce small, predictable changes in the output, unlike a calculation that can blow up from rounding) for short distances because it works with sin² of half-angles; the older spherical law of cosines, d = R arccos(...), takes arccos of a value extremely close to 1 for nearby points and loses most of its precision in double arithmetic. Clamping a to at most 1 guards Haversine against rounding just above 1 for nearly antipodal points.
• Vincenty convergence: 1 to 4 iterations here. Near-antipodal pairs can need many iterations or fail, which the code reports as NaN instead of looping forever.

When a planar approximation is acceptable in production

Use it when all four hold: the search radius is metropolitan (under about 50 km); the cos(latitude) factor is applied (precomputed once per city or per request); longitude differences are wrapped at the antimeridian; and the service does not operate near the poles (above about 80 degrees, switch to Haversine or a sphere-native index). Ride matching meets all four almost everywhere.

Effect on ranking, latency and system-level optimizations

• Ranking: the final ranking should use estimated travel time from a routing engine (a service, such as Google Maps or OSRM, that computes real driving routes and travel times over the actual road network); geometric distance only picks the candidate set (say the 50 nearest), so its tiny errors are irrelevant as long as it does not reorder candidates across the cut-off. The cos(latitude) factor is what protects the ordering.
• Latency: per pair, planar is a handful of multiplies, Haversine several transcendental functions (sin, cos, arcsin: functions with no finite formula of additions and multiplications, so a CPU evaluates them with a more expensive series or lookup), Vincenty an iterative loop with data-dependent branches. At 1,000 match requests per second x 500 candidates each, that is 500,000 distance evaluations per second: trivial with planar or vectorized Haversine, wasteful with Vincenty.
• Vectorization and batching: planar and Haversine vectorize cleanly (the same operations on every array element). Vincenty's variable iteration count breaks SIMD lanes, the parallel slots a single SIMD instruction advances in lockstep (all lanes wait for the slowest). Batch all candidates for one request into one array call.
• Approximations that preserve order: comparing squared planar distances avoids the square root; comparing the Haversine intermediate value a (which grows with distance) avoids the arcsin and square root. Apply the full formula only to the final few.

Pitfalls

• Passing degrees into trigonometric functions that expect radians, or swapping latitude and longitude.
• Using Euclidean distance on raw degrees, which stretches east-west distances.
• Believing Vincenty is "more correct" for matching: it is more exact geometrically, but the road network error is far larger. Even at the 1 km range in the table above, Haversine differs from Vincenty by only about 0.28 m, while a single one-way street or a block a car cannot cut through can add a couple of hundred metres to the real route, a gap that is easily thousands of times the 0.28 m the formula choice costs you.