Haversine vs Equirectangular Approximation Error Budgets

Haversine is the distance function everybody reaches for and is nine times more expensive than the approximation that is more accurate over the ranges a geofence actually spans. That sounds wrong until the two error sources are separated: haversine is exact on a sphere and the Earth is not one, so it carries a fixed 0.5% ellipsoidal error at every range; the equirectangular approximation is exact on the ellipsoid locally and its error grows with range. Below the crossover the cheap function wins on both axes. This page sits under geodesic vs planar distance for fence tests within Spatial Indexing for Real-Time Checks.

Concept and specification #

The equirectangular approximation projects a small neighbourhood onto a plane by scaling longitude differences by the cosine of latitude, then applies Pythagoras:

$$x=\Delta \lambda \cos !\left(\frac{{\varphi }_1+{\varphi }_2}{2}\right),y=\Delta \varphi ,d=R\sqrt{x^2+y^2}$$

Its error comes from treating the meridian convergence as constant over the separation, so it grows as roughly the cube of the angular separation — negligible at hundreds of metres, meaningful at tens of kilometres. Haversine solves the spherical triangle exactly:

$$d=2R\arcsin \sqrt{{\sin }^2!\frac{\Delta \varphi }{2}+\cos {\varphi }_1\cos {\varphi }_2{\sin }^2!\frac{\Delta \lambda }{2}}$$

so its error is entirely the sphere-versus-ellipsoid discrepancy: a proportional 0.3–0.5% depending on latitude and bearing, which grows linearly with distance and never goes away.

That difference in growth is the whole story. One error is proportional to distance and the other is super-linear, so they cross:

The crossover is far further out than intuition suggests — around 200 km rather than the few kilometres most people assume — because haversine’s proportional error is large. For a geofencing workload, where the great majority of distance computations are between a device and a fence within a few kilometres, the equirectangular approximation is more accurate at a ninth of the cost, and haversine’s role is confined to coarse long-range filters where neither error matters.

The important caveat is that these figures use the mean Earth radius. Using a latitude-appropriate radius removes most of haversine’s error and moves the crossover in; using the local radius of curvature in the equirectangular formula removes most of its systematic bias too. Both refinements are one multiplication, and both are usually skipped.

Step-by-step implementation #

1. Precompute the cosine per fence, not per call. The cos(φ) term is the approximation’s only transcendental operation and it is constant for a given fence’s neighbourhood.

from __future__ import annotations
import math

# Radius of curvature terms for WGS84, good to a few parts per million.
_A = 6_378_137.0
_E2 = 6.694379990e-3

def metres_per_degree(lat_deg: float) -> tuple[float, float]:
    """(metres per degree of latitude, metres per degree of longitude) here."""
    phi = math.radians(lat_deg)
    s = math.sin(phi)
    w = math.sqrt(1.0 - _E2 * s * s)
    m_lat = (_A * (1.0 - _E2) / (w ** 3)) * math.pi / 180.0
    m_lon = (_A / w) * math.cos(phi) * math.pi / 180.0
    return m_lat, m_lon

class LocalScale:
    """Per-fence scale factors, computed once at authoring time."""
    __slots__ = ("lat0", "lon0", "m_lat", "m_lon")

    def __init__(self, lat0: float, lon0: float) -> None:
        self.lat0, self.lon0 = lat0, lon0
        self.m_lat, self.m_lon = metres_per_degree(lat0)

    def distance_m(self, lat: float, lon: float) -> float:
        dy = (lat - self.lat0) * self.m_lat
        dx = (lon - self.lon0) * self.m_lon
        return math.hypot(dx, dy)

    def distance2_m2(self, lat: float, lon: float) -> float:
        """Squared distance — skip the sqrt when only a comparison is needed."""
        dy = (lat - self.lat0) * self.m_lat
        dx = (lon - self.lon0) * self.m_lon
        return dx * dx + dy * dy

2. Use the ellipsoidal radii of curvature, not a single mean radius. metres_per_degree above uses the meridional and normal radii, which removes the approximation’s systematic bias at essentially no extra cost, since the scale factors are computed once per fence.

3. Compare squared distances wherever possible. A radius test needs no square root — comparing $d^2$ against $r^2$ is the same predicate for 40% less work, and the sqrt is often the dominant instruction in the function.

4. Set the range bound from the error budget, not from a habit. A single scale pair is accurate to a metre out to roughly 8 km at temperate latitudes; beyond that, either recompute the scale at the midpoint or switch functions. The bound belongs in a constant with a comment naming the accuracy it buys.

5. Vectorise for batches. Both functions are pure arithmetic over arrays; NumPy evaluates the equirectangular form over 4,000 points in 21 µs against 1.9 ms for a Python loop, and the same batching applies to haversine at proportionally higher cost.

6. Never use either function across the antimeridian without normalising the longitude difference — a discontinuity that both formulas share and that has its own treatment in antimeridian and pole-crossing geofences.

Benchmark and verification #

Measured over 50 million distance computations against Karney’s geodesic as ground truth, points drawn from a real fleet’s distribution (median separation 1.2 km):

The second row is the recommendation and it is better than every haversine variant on every column: 9.5× cheaper than haversine and two orders of magnitude more accurate over this distribution. The comparison between rows one and two is the cheapest improvement available anywhere in this pipeline — using the proper radii of curvature instead of a mean radius costs 3 ns and reduces the median error by a factor of 78.

The max-error column is the one to watch when adopting this. 1.4 m occurs at the tail of the separation distribution — the handful of comparisons at tens of kilometres — and whether that matters depends on what the comparison decides. For a fence-proximity filter it is irrelevant; for a boundary test on a small fence at long range, which should not be happening anyway, it is a signal that the scale pair is stale for that query.

Verify by sampling: for 10,000 random point pairs from production traffic, compute the difference against a geodesic reference and export it as a histogram. Anything above a metre indicates a scale factor being used far from where it was computed.

Failure modes and edge cases #

The polar case in the fifth row is a genuine limit rather than a tunable. As $\cos \varphi \to 0$ the longitude scale collapses, so a fixed longitude difference maps to an ever-smaller ground distance and the approximation’s assumption of a locally flat, uniformly scaled plane fails. Above about 75° the honest answer is to use a geodesic or an azimuthal projection centred at the pole; for a fleet operating in Svalbard or northern Alaska that is not an edge case but the normal condition.

One efficiency note worth stating because it is easy to miss. Once a fence carries a LocalScale, the natural next step is to store the fence’s geometry pre-scaled into those same metres — at which point the containment test is ordinary planar geometry and the scale multiplication happens once per query rather than once per vertex. That is exactly the arrangement described in the parent topic, and the approximation on this page is what makes it cheap enough to justify.

Related #

• Geodesic vs Planar Distance for Geofence Tests — the parent topic and the per-fence projection this feeds.
• Local Projection Selection for Metric Buffers — when a full projection is needed instead of a scale pair.
• Speed Gating to Suppress Impossible Jumps — the highest-frequency consumer of this function.