Measuring the World with JavaScript
If I asked you to measure two points on a table, you might get out a ruler or a tape measure to tell me the distance between the two. However, what if I asked you to measure the points between two capital cities? Suddenly it becomes more complicated, not only just logistically but also in terms of assumptions we have to make about that measurement.
One big change, apart from just scale, is that we now have to account for the roundness of the world. Thankfully mathematicians who were far smarter than me have come up with formulas that can account for this in various ways. One such formula is the haversine formula. This is credited as far back as 1801 to Spanish astronomer and mathematician José de Mendoza y Ríos.
The Haversine formula is useful in a lot of spherical trigonometry problems, as it allows a person to determine what is called a great-circle distance between two given points on a sphere. In particular relevance to what we are looking at today, it is useful for determining the approximate distance between two latitude and longitudes on the globe.
When it comes to translating this to code Chris Veness, a UK based programmer has written a healthy amount of code in the space - a lot of it is featured in the popular geospatial library Turf.js. Here's a stand-alone adaptation of some of his code for the Haversine formula, written in TypeScript:
export function haversineDistance(
pointOne: { lng: number; lat: number },
pointTwo: { lng: number; lat: number }
) {
const toRadians = (latOrLng: number) => (latOrLng * Math.PI) / 180;
const phiOne = toRadians(pointOne.lat);
const lambdaOne = toRadians(pointOne.lng);
const phiTwo = toRadians(pointTwo.lat);
const lambdaTwo = toRadians(pointTwo.lng);
const deltaPhi = phiTwo - phiOne;
const deltalambda = lambdaTwo - lambdaOne;
const a =
Math.sin(deltaPhi / 2) * Math.sin(deltaPhi / 2) +
Math.cos(phiOne) *
Math.cos(phiTwo) *
Math.sin(deltalambda / 2) *
Math.sin(deltalambda / 2);
const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
const radius = 6371e3;
const distance = radius * c;
return distance;
}The Haversine formula suffers from one core problem, which is that it assumes that the world is round. Unfortunately for us, the world is not actually round, it's closer to what we'd determine as an oblate spheroid; a sphere that is flattened at the poles. As the Haversine makes this assumption, it cannot be guaranteed correct to provide results with better than 0.5% accuracy.
This is where another formula comes in, Vincenty's inverse formula. The formula comes from Thaddeus Vincenty, a Polish American geodesist who was with the U.S. Air Force. Vincenty's inverse formula works with the idea that the Earth is an oblate-spheroid. This approach is more complex and is also iteration based, unlike the haversine formula.
Again, I've taken some of Chris's code and made it standalone (you can just copy and paste this function). It works on the WGS84 ellipsoid, as this is probably the most common use case. Lets take a look at the code, again in TypeScript:
export function inverseVincentyDistance(
pointOne: { lng: number; lat: number },
pointTwo: { lng: number; lat: number }
): number {
const toRadians = (latOrLng: number) => (latOrLng * Math.PI) / 180;
const phiOne = toRadians(pointOne.lat);
const lambda1 = toRadians(pointOne.lng);
const phiTwo = toRadians(pointTwo.lat);
const lambda2 = toRadians(pointTwo.lng);
const wgs84ellipsoid = {
a: 6378137,
b: 6356752.314245,
f: 1 / 298.257223563,
};
const { a, b, f } = wgs84ellipsoid;
const L = lambda2 - lambda1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanphi.
const tanU1 = (1 - f) * Math.tan(phiOne),
cosU1 = 1 / Math.sqrt(1 + tanU1 * tanU1),
sinU1 = tanU1 * cosU1;
const tanU2 = (1 - f) * Math.tan(phiTwo),
cosU2 = 1 / Math.sqrt(1 + tanU2 * tanU2),
sinU2 = tanU2 * cosU2;
const antipodal =
Math.abs(L) > Math.PI / 2 || Math.abs(phiTwo - phiOne) > Math.PI / 2;
let lambda = L,
sinLambda = null,
cosLambda = null; // lambda = difference in longitude on an auxiliary sphere
let sigma = antipodal ? Math.PI : 0,
sinSigma = 0,
cosSigma = antipodal ? -1 : 1,
sinSqsigma = null; // sigma = angular distance P₁ P₂ on the sphere
let cos2sigmaM = 1; // sigmaM = angular distance on the sphere from the equator to the midpoint of the line
let sinalpha = null,
cosSqAlpha = 1; // alpha = azimuth of the geodesic at the equator
let C = null;
let lambdaʹ = null,
iterations = 0;
do {
sinLambda = Math.sin(lambda);
cosLambda = Math.cos(lambda);
sinSqsigma =
cosU2 * sinLambda * (cosU2 * sinLambda) +
(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) *
(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
if (Math.abs(sinSqsigma) < Number.EPSILON) {
break; // co-incident/antipodal points (falls back on lambda/sigma = L)
}
sinSigma = Math.sqrt(sinSqsigma);
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
sigma = Math.atan2(sinSigma, cosSigma);
sinalpha = (cosU1 * cosU2 * sinLambda) / sinSigma;
cosSqAlpha = 1 - sinalpha * sinalpha;
cos2sigmaM =
cosSqAlpha != 0 ? cosSigma - (2 * sinU1 * sinU2) / cosSqAlpha : 0; // on equatorial line cos²alpha = 0 (§6)
C = (f / 16) * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
lambdaʹ = lambda;
lambda =
L +
(1 - C) *
f *
sinalpha *
(sigma +
C *
sinSigma *
(cos2sigmaM + C * cosSigma * (-1 + 2 * cos2sigmaM * cos2sigmaM)));
const iterationCheck = antipodal
? Math.abs(lambda) - Math.PI
: Math.abs(lambda);
if (iterationCheck > Math.PI) {
throw new Error("lambda > Math.PI");
}
} while (Math.abs(lambda - lambdaʹ) > 1e-12 && ++iterations < 1000);
if (iterations >= 1000) {
throw new Error("Vincenty formula failed to converge");
}
const uSq = (cosSqAlpha * (a * a - b * b)) / (b * b);
const A = 1 + (uSq / 16384) * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
const B = (uSq / 1024) * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
const deltaSigma =
B *
sinSigma *
(cos2sigmaM +
(B / 4) *
(cosSigma * (-1 + 2 * cos2sigmaM * cos2sigmaM) -
(B / 6) *
cos2sigmaM *
(-3 + 4 * sinSigma * sinSigma) *
(-3 + 4 * cos2sigmaM * cos2sigmaM)));
const distance = b * A * (sigma - deltaSigma); // distance = length of the geodesic
return distance;
}I'll be honest, Vincenty's formula is a bit lost on me, but it promises to provide more accurate results! In fact, its use is very common in geodesy because the expected accuracy is to 0.5mm on the Earth ellipsoid. Pretty cool right? The one downside is that it is computationally slower due to the increased number of operations and the iterative approach.
Hopefully, you've enjoyed this short blog post about how to measure distances on the Earth with JavaScript (...okay TypeScript). Let me know if you would be interested in follow-ups on similar topics!
Published