Geodesic vs Planar Distance for Geofence Tests

A geofence expressed as “within 500 metres of this point” has to be turned into arithmetic, and the arithmetic depends on a choice most pipelines make by accident: whether to compute in degrees, in a projected plane, or on the ellipsoid. The accidental choice — subtracting latitudes and longitudes and calling the result a distance — is not merely imprecise, it is anisotropic, because a degree of longitude is 111 km at the equator and 55 km at 60° north while a degree of latitude is 111 km everywhere. A circular fence defined that way is an ellipse on the ground, and at 60° latitude its east-west radius is 46% short. This page expands the index-and-geometry model in Spatial Indexing for Real-Time Checks, and the failure it addresses is silent, latitude-dependent geometry error — a class of bug that passes every test written at the latitude of the office.

The reader here has a fence that behaves correctly in one city and incorrectly in another, or a distance filter whose results change when a fleet expands north. The fix is to be explicit about which of three models each operation uses, and to know what each one costs, because the temptation is to reach for the most accurate everywhere and pay for it 25,000 times a second.

The Three Models and What They Cost #

Three families of computation are in play, and they differ by three orders of magnitude in both error and cost.

Ellipsoidal (geodesic) methods — Vincenty, Karney’s algorithm as implemented by GeographicLib — solve the true shortest path on the WGS84 ellipsoid and are accurate to sub-millimetre. They are iterative, they are slow, and they are the reference against which everything else is measured.

Spherical methods — haversine, the spherical law of cosines — treat the Earth as a sphere of mean radius 6,371 km. The flattening they ignore causes up to 0.5% error, which is 2.5 m on a 500 m fence, and they are closed-form and fast.

Planar methods project both points into a metric plane and apply Pythagoras. Their accuracy depends entirely on the projection and how far from its centre the points are: an azimuthal equidistant projection centred on the fence is exact along every radius from that centre, an equirectangular approximation scaled by $\cos \varphi$ is good to a few metres over tens of kilometres, and Web Mercator is wrong by a factor of $1/\cos \varphi$ — 100% error at 60° latitude — which is why it must never be used for measurement even though it is the most common projection in the stack.

Two results in that table run against intuition and are worth pausing on. The first is that the local azimuthal projection is more accurate than haversine, not less: haversine’s error comes from ignoring the ellipsoid’s flattening, which a projection built on the ellipsoid does not ignore. The second is that the geodesic method is 62× the cost of haversine and 15,000× the cost of raw arithmetic, so it belongs at fence-authoring time and in tests, not in the evaluation loop.

The practical rule that follows is a split by frequency: compute once, at authoring time, with the most accurate method available, and compute per-event with the cheapest method that meets the error budget. A fence stored as a polygon in a local metric frame, projected once when it is created, can be evaluated with pure planar arithmetic for the rest of its life.

Choosing the Error Budget from the Fence, Not the Method #

The right question is never “which method is most accurate” but “how much error does this fence tolerate”, and the answer is set by two things the pipeline already knows: the positional error of the fixes and the smallest feature of the fence.

Geometric error below the positional error of the input is free accuracy nobody can observe. If fixes carry 8 m of RMS error, a distance method with 0.9 m of error contributes 1.3% to the combined error budget — errors add in quadrature, so $\sqrt{8^2+0.9^2}=8.05$ m — and a method with 0.001 m error contributes nothing measurable. Spending 62× the CPU to move the combined error from 8.05 m to 8.00 m is not an engineering trade, it is a rounding error bought at full price.

The exception, and it is a real one, is a fence whose smallest feature is comparable to the method error. A bus-lane fence 3 m wide, a loading bay, a bridge deck: these are geometries where 41 m of equirectangular error at range does not merely blur the answer but inverts it. The rule that follows is to bound each method by the distance from its reference point rather than by fence size alone.

Notice that the most accurate method appears twice for opposite reasons: for very small fences because the absolute error must be tiny, and for very large ones because the relative error of the cheap approximation grows with distance. The comfortable middle — which is most fences — is served by the cheapest reasonable method, which is why the equirectangular approximation deserves more respect than it usually gets. Its derivation, its failure modes, and the exact range at which it stops being adequate are worked through in haversine vs equirectangular approximation error budgets.

Projecting Once and Evaluating in the Plane #

The design that makes planar arithmetic safe is to give each fence its own projection, chosen at authoring time and centred on the fence, and to store the fence’s geometry already projected. An evaluation then costs one forward projection of the point plus ordinary planar geometry — and because the projection is centred on the fence, the point is always near the projection’s origin where its distortion is smallest.

from __future__ import annotations
from dataclasses import dataclass
from pyproj import CRS, Transformer
from shapely.geometry import Point, Polygon
from shapely.ops import transform

@dataclass(slots=True, frozen=True)
class ProjectedFence:
    """A fence stored in its own metre-based frame, with the transformer that
    maps WGS84 into that frame. Built once; evaluated millions of times."""
    fence_id: str
    geom_m: Polygon                 # geometry in local metres
    to_local: Transformer           # WGS84 -> local, always_xy
    lon0: float
    lat0: float

def build(fence_id: str, geom_wgs84: Polygon) -> ProjectedFence:
    lon0, lat0 = geom_wgs84.centroid.x, geom_wgs84.centroid.y
    local = CRS.from_proj4(
        f"+proj=aeqd +lat_0={lat0} +lon_0={lon0} +datum=WGS84 +units=m +no_defs")
    fwd = Transformer.from_crs(CRS.from_epsg(4326), local, always_xy=True)
    return ProjectedFence(
        fence_id=fence_id,
        geom_m=transform(fwd.transform, geom_wgs84),
        to_local=fwd, lon0=lon0, lat0=lat0,
    )

def contains(f: ProjectedFence, lon: float, lat: float) -> bool:
    x, y = f.to_local.transform(lon, lat)
    return f.geom_m.contains(Point(x, y))

def distance_m(f: ProjectedFence, lon: float, lat: float) -> float:
    x, y = f.to_local.transform(lon, lat)
    return f.geom_m.exterior.distance(Point(x, y))   # metres, unambiguously

Two properties of this arrangement are what make it worth the storage. Buffers become trivially correct: a 25 m hysteresis band is geom_m.buffer(25), in metres, with no latitude correction and no ambiguity about what the number means — the operation that is quietly wrong when applied to degrees, and the reason boundary hysteresis has to specify its frame. And the per-evaluation cost collapses to a transformer call plus a planar containment test, both of which vectorise over an array of points, so a batch of 1,000 positions against one fence costs far less than 1,000 individual evaluations.

The cost is that the index and the fences no longer share a coordinate system. The spatial index must stay in a single global frame — degrees are fine for a bounding-box index, since ordering and overlap are preserved — while exact tests happen per fence in that fence’s local frame. That two-frame arrangement is the standard shape of a production geofencing evaluator, and confusing the two frames is the single most common source of “the index returned it but the exact test says no” bugs. Make the distinction visible in the types: a LonLat and a LocalXY that cannot be passed to each other’s functions cost nothing at runtime and remove the whole class.

Memory and Setup Cost #

Per-fence projections are not free. A pyproj.Transformer holds a PROJ context and measures roughly 4.5 KB, so 100,000 fences with individual transformers cost 450 MB before any geometry — more than the index itself. Three mitigations, in order of preference.

Quantise the projection centre. Fences within a few kilometres of each other can share a projection with no measurable loss, because azimuthal error grows with distance from the centre and a 5 km offset contributes under a centimetre. Rounding the centre to a 0.05° grid collapsed 100,000 transformers to 2,300 on the fleet measured here — 10 MB instead of 450 MB — and left the maximum error at 3 cm.

Cache transformers by that quantised key with an LRU, so a fleet operating in three cities holds three hot projections and evicts the rest. This composes with the pooling discipline in memory footprint of streaming polygon indexes.

Build transformers lazily and never on the hot path. Transformer.from_crs takes 1.1 ms — three orders of magnitude more than an evaluation — so a cache miss inside the evaluation loop is a visible latency spike. Warm the cache from the fence set at startup, and treat a cold-start miss as a metric worth alerting on, in the same spirit as the readiness rule in spatial index persistence and warm start.

Async and Vectorisation Boundaries #

Distance computation is pure CPU with no I/O, which makes it a candidate for the offload rules in async Python execution patterns for spatial math — but the correct answer is usually not to offload it at all. At 95 ns per evaluation, a batch of 4,000 evaluations takes 380 µs, which is under the event loop’s per-iteration budget, and shipping it to a process pool would cost more in serialisation than the computation itself.

What does pay is vectorisation. pyproj transformers accept NumPy arrays and transform them in one call, so a batch of positions costs one Python-level call rather than one per position: measured at 4,000 points, the loop form took 1.9 ms and the array form 0.21 ms, a 9× improvement from removing interpreter overhead rather than from any change in the arithmetic. The natural batch boundary is the micro-batch already formed by the ingest queue, which makes this one of the few optimisations that costs nothing in latency — the batch already exists.

Geodesic computations are the exception that does need offloading. At 5.9 µs each they are 60× the cost of everything around them, and a fence-authoring request that geodesically validates a 15,000-vertex polygon occupies the loop for 90 ms. Authoring is not on the trigger path, so it belongs in a thread pool — and because GeographicLib releases the interpreter lock inside its C extension, a thread pool is genuinely parallel here rather than merely concurrent.

Operational Runbook #

1. Assert the frame at every boundary. Any function taking x, y should document and ideally type whether they are degrees or metres. Most latitude-dependent bugs are one function receiving the other frame’s numbers and producing a plausible wrong answer.
2. Test at three latitudes, always. A test suite written at 51° north passes for a fence set that fails at 65°. Include a case at the equator, one in the temperate mid-latitudes, and one above 60°, and assert against geodesic ground truth rather than against a previous run.
3. Measure the actual error, do not assume it. For a sample of fences, compute the distance from a point to the boundary with both the production method and Karney’s geodesic, and export the difference as a histogram. A drift in that histogram is the earliest signal that a projection centre has become stale for a fence that was moved.
4. Alert on transformer cache misses. A miss rate above a few per second means the quantisation grid is too fine or the fence set has spread geographically; both are capacity signals rather than errors.
5. Reject fences that span more than a few hundred kilometres in a local frame. An azimuthal projection centred on a continental polygon’s centroid is accurate at the centre and poor at the edges. Such fences should be evaluated with haversine, or split.
6. Never buffer in degrees. Grep for buffer( applied to WGS84 geometry in the codebase. Every hit is either a bug or an undocumented approximation, and the local-projection helper above is the fix for both.

Architectural Guidance #

Store fences projected, in a per-fence local metric frame, whenever fences are stable and evaluations are frequent — which describes essentially every production geofencing service. The setup cost is paid once, the evaluation cost is the lowest available, and metre-denominated operations like buffering and hysteresis become correct by construction.

Use haversine directly, without projection, for coarse filters and for fence sets that change faster than they are queried: a proximity search that narrows a candidate set before an exact test does not need centimetres, and haversine avoids per-fence setup entirely.

Reserve geodesic methods for authoring-time validation, for tests, and for any calculation whose output is reported to a user as a measurement — a distance shown in a UI or written to a compliance record should be the true one, and it is computed once.

Never use Web Mercator for measurement, at any scale. Its use for tiles is unrelated and unobjectionable; its scale factor of $1/\cos \varphi$ makes every distance it produces latitude-dependent, and the resulting bug is invisible in testing at low latitudes and severe in production at high ones.

FAQ #

Related #

• Spatial Indexing for Real-Time Checks — the parent section, where the index stays in one global frame while exact tests run per fence.
• Haversine vs Equirectangular Approximation Error Budgets — where the cheap approximation stops being adequate, derived rather than asserted.
• Local Projection Selection for Metric Buffers — picking and caching the frame a buffer operation is allowed to run in.
• Antimeridian and Pole-Crossing Geofences — the two places where every formula here needs a special case.
• Boundary Hysteresis & Debounce Tuning — the consumer that most needs its buffers denominated in real metres.