Advanced topics — Vincenty, ECEF, and affine transforms

Vincenty's formulae (ellipsoidal accuracy)

Haversine assumes a sphere; Vincenty works on the WGS84 ellipsoid and is accurate to ~0.5 mm. It's iterative and slow, but used for surveying, aviation, and any time you need sub-meter accuracy over long distances.

The Earth is not a perfect sphere — it's an oblate spheroid, with a radius at the equator (~6,378 km) that is 21 km larger than at the poles (~6,357 km). This difference is small (~0.3%) but accumulates: over a 10,000 km flight, Haversine and Vincenty can disagree by up to ~70 m. For most web apps this is irrelevant. For geodetic survey or precision GPS work, it's not.

The forward and inverse Vincenty formulas are too long to reproduce here, but every major language has a library. Use it when accuracy matters more than speed.

ECEF coordinates

Earth-Centered, Earth-Fixed is a 3D Cartesian system with origin at Earth's center, X through (0, 0), Y through (0, 90° E), Z through the North Pole.

$$X = (N + h) \cdot \cos\varphi \cdot \cos\lambda$$ $$Y = (N + h) \cdot \cos\varphi \cdot \sin\lambda$$ $$Z = (N \cdot (1 - e^2) + h) \cdot \sin\varphi$$

Where:

• $h$ is the height above the ellipsoid
• $N = a / \sqrt{1 - e^2 \sin^2\varphi}$ is the radius of curvature in the prime vertical
• $e^2$ is the squared eccentricity (≈ 0.00669437999 for WGS84)
• $a$ is the equatorial radius

ECEF is essential for GPS, sensor fusion, and 3D map engines (Cesium, Google Earth).

The radius of curvature $N$ is latitude-dependent — it varies from $a$ at the equator (6,378,137 m) to $a/(1-e^2)$ at the poles (~6,399,594 m). This is what distinguishes the ellipsoidal formulas from the spherical ones.

Affine transforms (2D map view)

When you pan and zoom a map, you're applying a 2D affine transform to world pixel coordinates:

$$\begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} = \begin{bmatrix} a & b & t_x \\ c & d & t_y \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$$

For pure pan (no rotation, no scale change): $a = d = 1$, $b = c = 0$, $t_x/t_y$ = offset. For zoom: $a = d = \text{scale}$, plus translation to keep the zoom anchor stationary.

If your map rotates by angle α: $a = \cos\alpha$, $b = -\sin\alpha$, $c = \sin\alpha$, $d = \cos\alpha$.