This script calculates great-circle distances between the two points – that is, the shortest distance over the earth’s surface – using the ‘Haversine’ formula.
It assumes a spherical earth, ignoring ellipsoidal effects – which is accurate enough^{*} for most purposes… – giving an ‘as-the-crow-flies’ distance between the two points (ignoring any hills!).
Enter the co-ordinates into the text boxes to try it out. It accepts a variety of formats:
Lat 1: Long 1:
Lat 2: Long 2:
And you can see it on a map (thanks to the nice guys at Google Maps)
Haversine formula:
R = earth’s radius (mean radius = 6,371km) Δlat = lat_{2}− lat_{1} Δlong = long_{2}− long_{1} a = sin²(Δlat/2) + cos(lat_{1}).cos(lat_{2}).sin²(Δlong/2) c = 2.atan2(√a, √(1−a)) d = R.c
(Note that angles need to be in radians to pass to trig functions).
var R = 6371; // km var dLat = (lat2-lat1).toRad(); var dLon = (lon2-lon1).toRad(); var a = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) * Math.sin(dLon/2) * Math.sin(dLon/2); var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); var d = R * c;
The Haversine formula ‘remains particularly well-conditioned for numerical computation even at small distances’ – unlike calculations based on the spherical law of cosines. (It was published by R W Sinnott in Sky and Telescope, 1984, though has been known about for much longer; the ‘half-versed-sine’ is (1-cosθ)/2, or sin²(θ/2) – don’t ask, I’m not a mathematician).
In fact, when Sinnott devised the Haversine formula, computational precision was limited. Nowadays, JavaScript (and most modern computers) use IEEE 754 64-bit floating-point numbers, which provide 15 significant figures of precision. With this precision, the simple spherical law of cosines formula gives well-conditioned results down to distances as small as around 1 metre. In view of this it is probably worth, in most situations, using either the simpler law of cosines or the more accurate ellipsoidal Vincenty formula in preference to Haversine! (See notes below on the limitations in accuracy of the spherical model).
var R = 6371; // km var d = Math.acos(Math.sin(lat1)*Math.sin(lat2) + Math.cos(lat1)*Math.cos(lat2) *Math.cos(lon2-lon1)) * R;
var y = Math.sin(dLon) * Math.cos(lat2); var x = Math.cos(lat1)*Math.sin(lat2) - Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon); var brng = Math.atan2(y, x).toBrng();
Since atan2 returns values in the range -π ... +π (that is, -180° ... +180°), to normalise the result to a compass bearing (in the range 0° ... 360°, with -ve values transformed into the range 180° ... 360°), convert to degrees and then use (θ+360) % 360, where % is modulo.
This is the initial bearing which if followed in a straight line along a great-circle arc (orthodrome) will take you from the start point to the end point; in general, the bearing you are following will have varied by the time you get to the end point (if you were to go from say 35°N,45°E (Baghdad) to 35°N,135°E (Osaka), you would start on a bearing of 60° and end up on a bearing of 120°!).
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
var Bx = Math.cos(lat2) * Math.cos(dLon); var By = Math.cos(lat2) * Math.sin(dLon); var lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2), Math.sqrt((Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By ) ); var lon3 = lon1.toRad() + Math.atan2(By, Math.cos(lat1) + Bx);
Just as the initial bearing may vary from the final bearing, the midpoint may not be located half-way between latitudes/longitudes; the midpoint between 35°N,45°E and 35°N,135°E is around 45°N,90°E.
This page is steadily growing! Given a start point, initial bearing, and distance, this will calculate the destination point and final bearing travelling along a (shortest distance) great circle arc.
var lat2 = Math.asin( Math.sin(lat1)*Math.cos(d/R) + Math.cos(lat1)*Math.sin(d/R)*Math.cos(brng) ); var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(d/R)*Math.cos(lat1), Math.cos(d/R)-Math.sin(lat1)*Math.sin(lat2));
Here’s a new one: I’ve sometimes been asked about distance of a point from a great-circle path (sometimes called cross track error).
var dXt = Math.asin(Math.sin(d13/R)*Math.sin(brng13-brng12)) * R;
Here, the great-circle path is identified by a start point and an end point – depending on what initial data you’re working from, you can use the formulae above to obtain the relevant distance and bearings. The sign of d_{xt} tells you which side of the path the third point is on.
The along-track distance, from the start point to the closest point on the path to the third point, is
var dAt = Math.acos(Math.cos(d13/R)/Math.cos(dXt/R)) * R;
Also new: ‘Clairaut’s formula’ will give you the maximum latitude of a great circle path, given a bearing and latitude on the great circle:
var latMax = Math.acos(Math.abs(Math.sin(brng)*Math.cos(lat)));
A ‘rhumb line’ (or loxodrome) is a path of constant bearing, which crosses all meridians at the same angle.
Sailors used to (and sometimes still) navigate along rhumb lines since it is easier to follow a constant compass bearing than to continually adjust the bearing as is needed to follow a great circle, though they are normally longer than great-circle (orthodrome) routes. Rhumb lines are straight lines on a Mercator Projection map.
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4)); var q = (Math.abs(dLat) > 1e-10) ? dLat/dPhi : Math.cos(lat1); // if dLon over 180° take shorter rhumb across 180° meridian: if (Math.abs(dLon) > Math.PI) { dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon); } var d = Math.sqrt(dLat*dLat + q*q*dLon*dLon) * R; var brng = Math.atan2(dLon, dPhi);
Given a start point and a distance d along constant bearing θ, this will calculate the destination point. If you maintain a constant bearing along a rhumb line, you will gradually spiral in towards one of the poles.
lat2 = lat1 + d*Math.cos(brng); var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4)); var q = (Math.abs(lat2-lat1) > 1e-10) ? (lat2-lat1)/dPhi : Math.cos(lat1); var dLon = d*Math.sin(brng)/q; // check for some daft bugger going past the pole, normalise latitude if so if (Math.abs(lat2) > Math.PI/2) lat2 = lat2>0 ? Math.PI-lat2 : -(Math.PI-lat2); lon2 = (lon1+dLon+Math.PI)%(2*Math.PI) - Math.PI;
If you use Ordnance Survey Grid References, I have implemented a script for converting between Lat/Long & OS Grid References.
No, I’ve not included decimal minutes: a decimal system is easy, a sexagesimal system has merits, but mixing the two is a complete sow’s ear. Switch off the option on your GPS!