Movable Type Scripts
Calculate distance, bearing and more between two Latitude/Longitude points
Notes:
- Accuracy: since the earth is not quite a sphere, there are small errors in using spherical geometry; the earth is actually roughly ellipsoidal (or more precisely, oblate spheroidal) with a radius varying between about 6,378km (equatorial) and 6,357km (polar), and local radius of curvature varying from 6,336km (equatorial meridian) to 6,399km (polar). This means that errors from assuming spherical geometry might be up to 0.55% crossing the equator, though generally below 0.3%, depending on latitude and direction of travel. An accuracy of better than 3m in 1km is good enough for me, but if you want greater accuracy, you could use the Vincenty formula for calculating geodesic distances on ellipsoids, which gives results accurate to within 1mm. (Out of sheer perversity – I’ve never needed such accuracy – I looked up this formula and discovered the JavaScript implementation was simpler than I expected).
- Trig functions take arguments in radians, so latitude, longitude, and bearings in degrees (either decimal or degrees/minutes/seconds) need to be converted to radians, rad = π.deg/180. When converting radians back to degrees (deg = 180.rad/π), West is negative if using signed decimal degrees. For bearings, values in the range -π to +π [-180° to +180°] need to be converted to 0 to +2π [0°–360°]; this can be done by (brng+2.π)%2.π [brng+360)%360] where % is the modulo operator. View page source to see JavaScript functions to handle these conversions.
- The atan2() function widely used here takes two arguments, atan2(y, x), and computes the arc tangent of the ratio y/x. It is more flexible than atan(y/x), since it handles x=0, and it also returns values in all 4 quadrants -π to +π (the atan function returns values in the range -π/2 to +π/2).
- If you implement any formula involving atan2 in Microsoft Excel, you will need to reverse the arguments, as Excel has them the opposite way around from JavaScript – conventional order is atan2(y, x), but Excel uses atan2(x, y). To use atan2 in a (VBA) macro, you can use WorksheetFunction.Atan2().
- All bearings are with respect to true north, 0°=N, 90°=E, etc; if you are working from a compass, magnetic north varies from true north in a complex way around the earth, and the difference has to be compensated for by variances indicated on local maps.
- I learned a lot from the US Census Bureau GIS FAQ which is no longer available, so I’ve made a copy.
- Thanks to Ed Williams’ Aviation Formulary for many of the formulae.
- For miles, divide km by 1.609344
- For nautical miles, divide km by 1.852
Principal JavaScript functions for distance calculation are shown below; use ‘View Source’ to see remaining functions for other formulae shown above and how they are used from HTML forms. These functions should be simple to translate into other languages if required.
Update for returning visitors, April 2007: As this page has grown, it’s become somewhat unweildy, so I have now revised and rationalised the scripts to made the code clearer and more re-usable. All LatLon methods now expect and return numeric degrees; converting between radians, degrees & deg/min/sec I have made extensions of Number or String objects. The LatLon object is just a handy way of holding a (latititude,longitude) pair. If you can’t or don’t want to use JavaScript objects, functions could perhaps return a point as an array (e.g. return [lat, lon];).
You are welcome to re-use these scripts [without any warranty express or implied] provided you retain my copyright notice and when possible a link to my website (under a LGPL license). If you have any queries or find any problems, please contact me.
© 2002-2008 Chris Veness
/*
* Use Haversine formula to Calculate distance (in km) between two points specified by
* latitude/longitude (in numeric degrees)
*
* example usage from form:
* result.value = LatLon.distHaversine(lat1.value.parseDeg(), long1.value.parseDeg(),
* lat2.value.parseDeg(), long2.value.parseDeg());
* where lat1, long1, lat2, long2, and result are form fields
*/
LatLon.distHaversine = function(lat1, lon1, lat2, lon2) {
var R = 6371; // earth's mean radius in km
var dLat = (lat2-lat1).toRad();
var dLon = (lon2-lon1).toRad();
lat1 = lat1.toRad(), lat2 = lat2.toRad();
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(lat1) * Math.cos(lat2) *
Math.sin(dLon/2) * Math.sin(dLon/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;
return d;
}
/*
* ditto using Law of Cosines
*/
LatLon.distCosineLaw = function(lat1, lon1, lat2, lon2) {
var R = 6371; // earth's mean radius in km
var d = Math.acos(Math.sin(lat1.toRad())*Math.sin(lat2.toRad()) +
Math.cos(lat1.toRad())*Math.cos(lat2.toRad())*Math.cos((lon2-lon1).toRad())) * R;
return d;
}
/*
* calculate (initial) bearing between two points
*
* from: Ed Williams' Aviation Formulary, http://williams.best.vwh.net/avform.htm#Crs
*/
LatLon.bearing = function(lat1, lon1, lat2, lon2) {
lat1 = lat1.toRad(); lat2 = lat2.toRad();
var dLon = (lon2-lon1).toRad();
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);
return Math.atan2(y, x).toBrng();
}
/*
* calculate destination point given start point, initial bearing (deg) and distance (km)
* see http://williams.best.vwh.net/avform.htm#LL
*/
LatLon.prototype.destPoint = function(brng, d) {
var R = 6371; // earth's mean radius in km
var lat1 = this.lat.toRad(), lon1 = this.lon.toRad();
brng = brng.toRad();
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));
lon2 = (lon2+Math.PI)%(2*Math.PI) - Math.PI; // normalise to -180...+180
if (isNaN(lat2) || isNaN(lon2)) return null;
return new LatLon(lat2.toDeg(), lon2.toDeg());
}
/*
* construct a LatLon object: arguments in numeric degrees
*
* note all LatLong methods expect & return numeric degrees (for lat/long & for bearings)
*/
function LatLon(lat, lon) {
this.lat = lat;
this.lon = lon;
}
/*
* represent point {lat, lon} in standard representation
*/
LatLon.prototype.toString = function() {
return this.lat.toLat() + ', ' + this.lon.toLon();
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
// extend String object with method for parsing degrees or lat/long values to numeric degrees
//
// this is very flexible on formats, allowing signed decimal degrees, or deg-min-sec suffixed by
// compass direction (NSEW). A variety of separators are accepted (eg 3º 37' 09"W) or fixed-width
// format without separators (eg 0033709W). Seconds and minutes may be omitted. (Minimal validation
// is done).
String.prototype.parseDeg = function() {
if (!isNaN(this)) return Number(this); // signed decimal degrees without NSEW
var degLL = this.replace(/^-/,'').replace(/[NSEW]/i,''); // strip off any sign or compass dir'n
var dms = degLL.split(/[^0-9.]+/); // split out separate d/m/s
for (var i in dms) if (dms[i]=='') dms.splice(i,1); // remove empty elements (see note below)
switch (dms.length) { // convert to decimal degrees...
case 3: // interpret 3-part result as d/m/s
var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600; break;
case 2: // interpret 2-part result as d/m
var deg = dms[0]/1 + dms[1]/60; break;
case 1: // decimal or non-separated dddmmss
if (/[NS]/i.test(this)) degLL = '0' + degLL; // - normalise N/S to 3-digit degrees
var deg = dms[0].slice(0,3)/1 + dms[0].slice(3,5)/60 + dms[0].slice(5)/3600; break;
default: return NaN;
}
if (/^-/.test(this) || /[WS]/i.test(this)) deg = -deg; // take '-', west and south as -ve
return deg;
}
// note: whitespace at start/end will split() into empty elements (except in IE)
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
// extend Number object with methods for converting degrees/radians
Number.prototype.toRad = function() { // convert degrees to radians
return this * Math.PI / 180;
}
Number.prototype.toDeg = function() { // convert radians to degrees (signed)
return this * 180 / Math.PI;
}
Number.prototype.toBrng = function() { // convert radians to degrees (as bearing: 0...360)
return (this.toDeg()+360) % 360;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
// extend Number object with methods for presenting bearings & lat/longs
Number.prototype.toDMS = function() { // convert numeric degrees to deg/min/sec
var d = Math.abs(this); // (unsigned result ready for appending compass dir'n)
d += 1/7200; // add ½ second for rounding
var deg = Math.floor(d);
var min = Math.floor((d-deg)*60);
var sec = Math.floor((d-deg-min/60)*3600);
// add leading zeros if required
if (deg<100) deg = '0' + deg; if (deg<10) deg = '0' + deg;
if (min<10) min = '0' + min;
if (sec<10) sec = '0' + sec;
return deg + '\u00B0' + min + '\u2032' + sec + '\u2033';
}
Number.prototype.toLat = function() { // convert numeric degrees to deg/min/sec latitude
return this.toDMS().slice(1) + (this<0 ? 'S' : 'N'); // knock off initial '0' for lat!
}
Number.prototype.toLon = function() { // convert numeric degrees to deg/min/sec longitude
return this.toDMS() + (this>0 ? 'E' : 'W');
}
Number.prototype.toPrecision = function(fig) { // override toPrecision method with one which displays
if (this == 0) return 0; // trailing zeros in place of exponential notation
var scale = Math.ceil(Math.log(this)*Math.LOG10E);
var mult = Math.pow(10, fig-scale);
return Math.round(this*mult)/mult;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */