mt-latlon.js

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/*  Latitude/longitude spherical geodesy formulae & scripts (c) Chris Veness 2002-2012            */
/*   - www.movable-type.co.uk/scripts/latlong.html                                                */
/*                                                                                                */
/*  Sample usage:                                                                                 */
/*    var p1 = new LatLon(51.5136, -0.0983);                                                      */
/*    var p2 = new LatLon(51.4778, -0.0015);                                                      */
/*    var dist = p1.distanceTo(p2);          // in km                                             */
/*    var brng = p1.bearingTo(p2);           // in degrees clockwise from north                   */
/*    ... etc                                                                                     */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/*  Note that minimal error checking is performed in this example code!                           */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */


(function (root) {
  'use strict';


  var wrapper = function(Geo) {

    /**
     * Creates a point on the earth's surface at the supplied latitude / longitude
     *
     * @constructor
     * @param {Number} lat: latitude in numeric degrees
     * @param {Number} lon: longitude in numeric degrees
     * @param {Number} [rad=6371]: radius of earth if different value is required from standard 6,371km
     */
    function LatLon(lat, lon, rad) {
      if (typeof(rad) == 'undefined') rad = 6371;  // earth's mean radius in km
      // only accept numbers or valid numeric strings
      this._lat = typeof(lat)=='number' ? lat : typeof(lat)=='string' && lat.trim()!='' ? +lat : NaN;
      this._lon = typeof(lon)=='number' ? lon : typeof(lon)=='string' && lon.trim()!='' ? +lon : NaN;
      this._radius = typeof(rad)=='number' ? rad : typeof(rad)=='string' && trim(lon)!='' ? +rad : NaN;
    }


    /**
     * Returns the distance from this point to the supplied point, in km
     * (using Haversine formula)
     *
     * from: Haversine formula - R. W. Sinnott, "Virtues of the Haversine",
     *       Sky and Telescope, vol 68, no 2, 1984
     *
     * @param   {LatLon} point: Latitude/longitude of destination point
     * @param   {Number} [precision=4]: no of significant digits to use for returned value
     * @returns {Number} Distance in km between this point and destination point
     */
    LatLon.prototype.distanceTo = function(point, precision) {
      // default 4 sig figs reflects typical 0.3% accuracy of spherical model
      if (typeof precision == 'undefined') precision = 4;

      var R = this._radius;
      var lat1 = toRad(this._lat), lon1 = toRad(this._lon);
      var lat2 = toRad(point._lat), lon2 = toRad(point._lon);
      var dLat = lat2 - lat1;
      var dLon = lon2 - lon1;

      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 toPrecisionFixed(d, precision);
    }


    /**
     * Returns the (initial) bearing from this point to the supplied point, in degrees
     *   see http://williams.best.vwh.net/avform.htm#Crs
     *
     * @param   {LatLon} point: Latitude/longitude of destination point
     * @returns {Number} Initial bearing in degrees from North
     */
    LatLon.prototype.bearingTo = function(point) {
      var lat1 = toRad(this._lat), lat2 = toRad(point._lat);
      var dLon = toRad(point._lon-this._lon);

      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);

      return (toDeg(brng)+360) % 360;
    }


    /**
     * Returns final bearing arriving at supplied destination point from this point; the final bearing
     * will differ from the initial bearing by varying degrees according to distance and latitude
     *
     * @param   {LatLon} point: Latitude/longitude of destination point
     * @returns {Number} Final bearing in degrees from North
     */
    LatLon.prototype.finalBearingTo = function(point) {
      // get initial bearing from supplied point back to this point...
      var lat1 = toRad(point._lat), lat2 = toRad(this._lat);
      var dLon = toRad(this._lon-point._lon);

      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);

      // ... & reverse it by adding 180°
      return (toDeg(brng)+180) % 360;
    }


    /**
     * Returns the midpoint between this point and the supplied point.
     *   see http://mathforum.org/library/drmath/view/51822.html for derivation
     *
     * @param   {LatLon} point: Latitude/longitude of destination point
     * @returns {LatLon} Midpoint between this point and the supplied point
     */
    LatLon.prototype.midpointTo = function(point) {
      var lat1 = toRad(this._lat), lon1 = toRad(this._lon);
      var lat2 = toRad(point._lat);
      var dLon = toRad(point._lon-this._lon);

      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 + Math.atan2(By, Math.cos(lat1) + Bx);
      var lon3 = (lon3+3*Math.PI) % (2*Math.PI) - Math.PI;  // normalise to -180..+180º

      return toLatLonDeg(lat3, lon3, this._radius);
    }


    /**
     * Returns the destination point from this point having travelled the given distance (in km) on the
     * given initial bearing (bearing may vary before destination is reached)
     *
     *   see http://williams.best.vwh.net/avform.htm#LL
     *
     * @param   {Number} brng: Initial bearing in degrees
     * @param   {Number} dist: Distance in km
     * @returns {LatLon} Destination point
     */
    LatLon.prototype.destinationPoint = function(brng, dist) {
      dist = typeof(dist)=='number' ? dist : typeof(dist)=='string' && dist.trim()!='' ? +dist : NaN;
      dist = dist/this._radius;  // convert dist to angular distance in radians
      brng = toRad(brng);  //
      var lat1 = toRad(this._lat), lon1 = toRad(this._lon);

      var lat2 = Math.asin( Math.sin(lat1)*Math.cos(dist) +
                            Math.cos(lat1)*Math.sin(dist)*Math.cos(brng) );
      var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(dist)*Math.cos(lat1),
                                   Math.cos(dist)-Math.sin(lat1)*Math.sin(lat2));
      lon2 = (lon2+3*Math.PI) % (2*Math.PI) - Math.PI;  // normalise to -180..+180º

      return toLatLonDeg(lat2, lon2, this._radius);
    }


    /**
     * Returns the point of intersection of two paths defined by point and bearing
     *
     *   see http://williams.best.vwh.net/avform.htm#Intersection
     *
     * @param   {LatLon} p1: First point
     * @param   {Number} brng1: Initial bearing from first point
     * @param   {LatLon} p2: Second point
     * @param   {Number} brng2: Initial bearing from second point
     * @returns {LatLon} Destination point (null if no unique intersection defined)
     */
    LatLon.intersection = function(p1, brng1, p2, brng2) {
      brng1 = typeof brng1 == 'number' ? brng1 : typeof brng1 == 'string' && trim(brng1)!='' ? +brng1 : NaN;
      brng2 = typeof brng2 == 'number' ? brng2 : typeof brng2 == 'string' && trim(brng2)!='' ? +brng2 : NaN;
      var lat1 = toRad(p1._lat), lon1 = toRad(p1._lon);
      var lat2 = toRad(p2._lat), lon2 = toRad(p2._lon);
      var brng13 = toRad(brng1), brng23 = toRad(brng2);
      var dLat = lat2-lat1, dLon = lon2-lon1;

      var dist12 = 2*Math.asin( Math.sqrt( Math.sin(dLat/2)*Math.sin(dLat/2) +
        Math.cos(lat1)*Math.cos(lat2)*Math.sin(dLon/2)*Math.sin(dLon/2) ) );
      if (dist12 == 0) return null;

      // initial/final bearings between points
      var brngA = Math.acos( ( Math.sin(lat2) - Math.sin(lat1)*Math.cos(dist12) ) /
        ( Math.sin(dist12)*Math.cos(lat1) ) );
      if (isNaN(brngA)) brngA = 0;  // protect against rounding
      var brngB = Math.acos( ( Math.sin(lat1) - Math.sin(lat2)*Math.cos(dist12) ) /
        ( Math.sin(dist12)*Math.cos(lat2) ) );

      var brng12, brng21;
      if (Math.sin(lon2-lon1) > 0) {
        brng12 = brngA;
        brng21 = 2*Math.PI - brngB;
      } else {
        brng12 = 2*Math.PI - brngA;
        brng21 = brngB;
      }

      var alpha1 = (brng13 - brng12 + Math.PI) % (2*Math.PI) - Math.PI;  // angle 2-1-3
      var alpha2 = (brng21 - brng23 + Math.PI) % (2*Math.PI) - Math.PI;  // angle 1-2-3

      if (Math.sin(alpha1)==0 && Math.sin(alpha2)==0) return null;  // infinite intersections
      if (Math.sin(alpha1)*Math.sin(alpha2) < 0) return null;       // ambiguous intersection

      //alpha1 = Math.abs(alpha1);
      //alpha2 = Math.abs(alpha2);
      // ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?

      var alpha3 = Math.acos( -Math.cos(alpha1)*Math.cos(alpha2) +
                           Math.sin(alpha1)*Math.sin(alpha2)*Math.cos(dist12) );
      var dist13 = Math.atan2( Math.sin(dist12)*Math.sin(alpha1)*Math.sin(alpha2),
                           Math.cos(alpha2)+Math.cos(alpha1)*Math.cos(alpha3) )
      var lat3 = Math.asin( Math.sin(lat1)*Math.cos(dist13) +
                        Math.cos(lat1)*Math.sin(dist13)*Math.cos(brng13) );
      var dLon13 = Math.atan2( Math.sin(brng13)*Math.sin(dist13)*Math.cos(lat1),
                           Math.cos(dist13)-Math.sin(lat1)*Math.sin(lat3) );
      var lon3 = lon1+dLon13;
      lon3 = (lon3+3*Math.PI) % (2*Math.PI) - Math.PI;  // normalise to -180..+180º

      return toLatLonDeg(lat3, lon3, this._radius);
    }


    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

    /**
     * Returns the distance from this point to the supplied point, in km, travelling along a rhumb line
     *
     *   see http://williams.best.vwh.net/avform.htm#Rhumb
     *
     * @param   {LatLon} point: Latitude/longitude of destination point
     * @returns {Number} Distance in km between this point and destination point
     */
    LatLon.prototype.rhumbDistanceTo = function(point) {
      var R = this._radius;
      var lat1 = toRad(this._lat), lat2 = toRad(point._lat);
      var dLat = toRad(point._lat-this._lat);
      var dLon = toRad(Math.abs(point._lon-this._lon));

      var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
      var q = (isFinite(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1);  // E-W line gives dPhi=0

      // if dLon over 180° take shorter rhumb across anti-meridian:
      if (Math.abs(dLon) > Math.PI) {
        dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);
      }

      var dist = Math.sqrt(dLat*dLat + q*q*dLon*dLon) * R;

      return toPrecisionFixed(dist, 4);  // 4 sig figs reflects typical 0.3% accuracy of spherical model
    }

    /**
     * Returns the bearing from this point to the supplied point along a rhumb line, in degrees
     *
     * @param   {LatLon} point: Latitude/longitude of destination point
     * @returns {Number} Bearing in degrees from North
     */
    LatLon.prototype.rhumbBearingTo = function(point) {
      var lat1 = toRad(this._lat), lat2 = toRad(point._lat);
      var dLon = toRad(point._lon-this._lon);

      var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
      if (Math.abs(dLon) > Math.PI) dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);
      var brng = Math.atan2(dLon, dPhi);

      return (toDeg(brng)+360) % 360;
    }

    /**
     * Returns the destination point from this point having travelled the given distance (in km) on the
     * given bearing along a rhumb line
     *
     * @param   {Number} brng: Bearing in degrees from North
     * @param   {Number} dist: Distance in km
     * @returns {LatLon} Destination point
     */
    LatLon.prototype.rhumbDestinationPoint = function(brng, dist) {
      var R = this._radius;
      var d = parseFloat(dist)/R;  // d = angular distance covered on earth’s surface
      var lat1 = toRad(this._lat), lon1 = toRad(this._lon);
      brng = toRad(brng);

      var dLat = d*Math.cos(brng);
      // nasty kludge to overcome ill-conditioned results around parallels of latitude:
      if (Math.abs(dLat) < 1e-10) dLat = 0; // dLat < 1 mm

      var lat2 = lat1 + dLat;
      var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
      var q = (isFinite(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1);  // E-W line gives dPhi=0
      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;

      var lon2 = (lon1+dLon+3*Math.PI)%(2*Math.PI) - Math.PI;

      return toLatLonDeg(lat2, lon2, this._radius);
    }

    /**
     * Returns the loxodromic midpoint (along a rhumb line) between this point and the supplied point.
     *   see http://mathforum.org/kb/message.jspa?messageID=148837
     *
     * @param   {LatLon} point: Latitude/longitude of destination point
     * @returns {LatLon} Midpoint between this point and the supplied point
     */
    LatLon.prototype.rhumbMidpointTo = function(point) {
      var lat1 = toRad(this._lat), lon1 = toRad(this._lon);
      var lat2 = toRad(point._lat), lon2 = toRad(point._lon);

      if (Math.abs(lon2-lon1) > Math.PI) lon1 += 2*Math.PI; // crossing anti-meridian

      var lat3 = (lat1+lat2)/2;
      var f1 = Math.tan(Math.PI/4 + lat1/2);
      var f2 = Math.tan(Math.PI/4 + lat2/2);
      var f3 = Math.tan(Math.PI/4 + lat3/2);
      var lon3 = ( (lon2-lon1)*Math.log(f3) + lon1*Math.log(f2) - lon2*Math.log(f1) ) / Math.log(f2/f1);

      if (!isFinite(lon3)) lon3 = (lon1+lon2)/2; // parallel of latitude

      var lon3 = (lon3+3*Math.PI) % (2*Math.PI) - Math.PI;  // normalise to -180..+180º

      return toLatLonDeg(lat3, lon3, this._radius);
    }


    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */


    /**
     * Returns the latitude of this point; signed numeric degrees if no format, otherwise format & dp
     * as per Geo.toLat()
     *
     * @param   {String} [format]: Return value as 'd', 'dm', 'dms'
     * @param   {Number} [dp=0|2|4]: No of decimal places to display
     * @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
     */
    LatLon.prototype.lat = function(format, dp) {
      if (typeof format == 'undefined') return this._lat;

      return Geo.toLat(this._lat, format, dp);
    }

    /**
     * Returns the longitude of this point; signed numeric degrees if no format, otherwise format & dp
     * as per Geo.toLon()
     *
     * @param   {String} [format]: Return value as 'd', 'dm', 'dms'
     * @param   {Number} [dp=0|2|4]: No of decimal places to display
     * @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
     */
    LatLon.prototype.lon = function(format, dp) {
      if (typeof format == 'undefined') return this._lon;

      return Geo.toLon(this._lon, format, dp);
    }

    /**
     * Returns a string representation of this point; format and dp as per lat()/lon()
     *
     * @param   {String} [format]: Return value as 'd', 'dm', 'dms'
     * @param   {Number} [dp=0|2|4]: No of decimal places to display
     * @returns {String} Comma-separated latitude/longitude
     */
    LatLon.prototype.toString = function(format, dp) {
      if (typeof format == 'undefined') format = 'dms';

      return Geo.toLat(this._lat, format, dp) + ', ' + Geo.toLon(this._lon, format, dp);
    }

    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

    // ---- helper functions for converting degrees/radians

    /** Converts numeric degrees to radians */
    function toRad(value) {
      return value * Math.PI / 180;
    }

    /** Converts radians to numeric (signed) degrees */
    function toDeg(value) {
      return value * 180 / Math.PI;
    }

    /**
     * Formats the significant digits of a number, using only fixed-point notation (no exponential)
     *
     * @param   {Number} value: The number to convert
     * @param   {Number} precision: Number of significant digits to appear in the returned string
     * @returns {String} A string representation of number which contains precision significant digits
     */
    function toPrecisionFixed(value, precision) {
      // use standard toPrecision method
      var n = value.toPrecision(precision);

      // ... but replace +ve exponential format with trailing zeros
      n = n.replace(/(.+)e\+(.+)/, function(n, sig, exp) {
        sig = sig.replace(/\./, '');       // remove decimal from significand
        var l = sig.length - 1;
        while (exp-- > l) sig = sig + '0'; // append zeros from exponent
        return sig;
      });

      // ... and replace -ve exponential format with leading zeros
      n = n.replace(/(.+)e-(.+)/, function(n, sig, exp) {
        sig = sig.replace(/\./, '');       // remove decimal from significand
        while (exp-- > 1) sig = '0' + sig; // prepend zeros from exponent
        return '0.' + sig;
      });

      return n;
    }

    function toLatLonDeg(latRad, lonRad, radius) {
      return new LatLon(toDeg(latRad), toDeg(lonRad), radius)
    }

    return LatLon;
  };

  // Add support for AMD, Node and plain JS
  if (typeof exports === 'object') {
      module.exports = wrapper;
  } else if (typeof define === 'function' && define.amd) {
      define(wrapper);
  } else {
      root.LatLon = wrapper(Geo);
  }

})(this);