How to calculate the latlng of a point a certain distance away from another?

Solution 1:

We will need a method that returns the destination point when given a bearing and the distance travelled from a source point. Luckily, there is a very good JavaScript implementation by Chris Veness at Calculate distance, bearing and more between Latitude/Longitude points.

The following has been adapted to work with the google.maps.LatLng class:

Number.prototype.toRad = function() {
   return this * Math.PI / 180;
}

Number.prototype.toDeg = function() {
   return this * 180 / Math.PI;
}

google.maps.LatLng.prototype.destinationPoint = function(brng, dist) {
   dist = dist / 6371;  
   brng = brng.toRad();  

   var lat1 = this.lat().toRad(), lon1 = this.lng().toRad();

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

   if (isNaN(lat2) || isNaN(lon2)) return null;

   return new google.maps.LatLng(lat2.toDeg(), lon2.toDeg());
}

You would simply use it as follows:

var pointA = new google.maps.LatLng(25.48, -71.26); 
var radiusInKm = 10;

var pointB = pointA.destinationPoint(90, radiusInKm);

Here is a complete example using Google Maps API v3:

<!DOCTYPE html>
<html> 
<head> 
   <meta http-equiv="content-type" content="text/html; charset=UTF-8"/> 
   <title>Google Maps Geometry</title> 
   <script src="http://maps.google.com/maps/api/js?sensor=false" 
           type="text/javascript"></script> 
</head> 
<body> 
   <div id="map" style="width: 400px; height: 300px"></div> 

   <script type="text/javascript"> 
      Number.prototype.toRad = function() {
         return this * Math.PI / 180;
      }

      Number.prototype.toDeg = function() {
         return this * 180 / Math.PI;
      }

      google.maps.LatLng.prototype.destinationPoint = function(brng, dist) {
         dist = dist / 6371;  
         brng = brng.toRad();  

         var lat1 = this.lat().toRad(), lon1 = this.lng().toRad();

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

         if (isNaN(lat2) || isNaN(lon2)) return null;

         return new google.maps.LatLng(lat2.toDeg(), lon2.toDeg());
      }

      var pointA = new google.maps.LatLng(40.70, -74.00);   // Circle center
      var radius = 10;                                      // 10km

      var mapOpt = { 
         mapTypeId: google.maps.MapTypeId.TERRAIN,
         center: pointA,
         zoom: 10
      };

      var map = new google.maps.Map(document.getElementById("map"), mapOpt);

      // Draw the circle
      new google.maps.Circle({
         center: pointA,
         radius: radius * 1000,       // Convert to meters
         fillColor: '#FF0000',
         fillOpacity: 0.2,
         map: map
      });

      // Show marker at circle center
      new google.maps.Marker({
         position: pointA,
         map: map
      });

      // Show marker at destination point
      new google.maps.Marker({
         position: pointA.destinationPoint(90, radius),
         map: map
      });
   </script> 
</body> 
</html>

Screenshot:

Google Maps Geometry

UPDATE:

In reply to Paul's comment below, this is what happens when the circle wraps around one of the poles.

Plotting pointA near the north pole, with a radius of 1,000km:

  var pointA = new google.maps.LatLng(85, 0);   // Close to north pole
  var radius = 1000;                            // 1000km

Screenshot for pointA.destinationPoint(90, radius):

Close to north pole

Solution 2:

To calculate a lat,long point at a given bearing and distance from another you can use google´s JavaScript implementation:

var pointA = new google.maps.LatLng(25.48, -71.26); 
var distance = 10; // 10 metres
var bearing = 90; // 90 degrees
var pointB = google.maps.geometry.spherical.computeOffset(pointA, distance, bearing);

See https://developers.google.com/maps/documentation/javascript/reference#spherical For documentation

Solution 3:

If you are after the distance between 2 lat/lng points across the earths surface then you can find the javascript here:

http://www.movable-type.co.uk/scripts/latlong-vincenty.html

This is the same formula used in android API at android.location.Location::distanceTo

You can easily convert the code from javascript to java.

If you want to calculate the destination point given the start point, bearing and distance, then you need this method:

http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html

Here are the formulae in java:

public class LatLngUtils {

  /**
   * @param lat1
   *          Initial latitude
   * @param lon1
   *          Initial longitude
   * @param lat2
   *          destination latitude
   * @param lon2
   *          destination longitude
   * @param results
   *          To be populated with the distance, initial bearing and final
   *          bearing
   */

  public static void computeDistanceAndBearing(double lat1, double lon1,
      double lat2, double lon2, double results[]) {
    // Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
    // using the "Inverse Formula" (section 4)

    int MAXITERS = 20;
    // Convert lat/long to radians
    lat1 *= Math.PI / 180.0;
    lat2 *= Math.PI / 180.0;
    lon1 *= Math.PI / 180.0;
    lon2 *= Math.PI / 180.0;

    double a = 6378137.0; // WGS84 major axis
    double b = 6356752.3142; // WGS84 semi-major axis
    double f = (a - b) / a;
    double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);

    double L = lon2 - lon1;
    double A = 0.0;
    double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
    double U2 = Math.atan((1.0 - f) * Math.tan(lat2));

    double cosU1 = Math.cos(U1);
    double cosU2 = Math.cos(U2);
    double sinU1 = Math.sin(U1);
    double sinU2 = Math.sin(U2);
    double cosU1cosU2 = cosU1 * cosU2;
    double sinU1sinU2 = sinU1 * sinU2;

    double sigma = 0.0;
    double deltaSigma = 0.0;
    double cosSqAlpha = 0.0;
    double cos2SM = 0.0;
    double cosSigma = 0.0;
    double sinSigma = 0.0;
    double cosLambda = 0.0;
    double sinLambda = 0.0;

    double lambda = L; // initial guess
    for (int iter = 0; iter < MAXITERS; iter++) {
      double lambdaOrig = lambda;
      cosLambda = Math.cos(lambda);
      sinLambda = Math.sin(lambda);
      double t1 = cosU2 * sinLambda;
      double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
      double sinSqSigma = t1 * t1 + t2 * t2; // (14)
      sinSigma = Math.sqrt(sinSqSigma);
      cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
      sigma = Math.atan2(sinSigma, cosSigma); // (16)
      double sinAlpha = (sinSigma == 0) ? 0.0 : cosU1cosU2 * sinLambda
          / sinSigma; // (17)
      cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
      cos2SM = (cosSqAlpha == 0) ? 0.0 : cosSigma - 2.0 * sinU1sinU2
          / cosSqAlpha; // (18)

      double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
      A = 1 + (uSquared / 16384.0) * // (3)
          (4096.0 + uSquared * (-768 + uSquared * (320.0 - 175.0 * uSquared)));
      double B = (uSquared / 1024.0) * // (4)
          (256.0 + uSquared * (-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
      double C = (f / 16.0) * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
      double cos2SMSq = cos2SM * cos2SM;
      deltaSigma = B
          * sinSigma
          * // (6)
          (cos2SM + (B / 4.0)
              * (cosSigma * (-1.0 + 2.0 * cos2SMSq) - (B / 6.0) * cos2SM
                  * (-3.0 + 4.0 * sinSigma * sinSigma)
                  * (-3.0 + 4.0 * cos2SMSq)));

      lambda = L
          + (1.0 - C)
          * f
          * sinAlpha
          * (sigma + C * sinSigma
              * (cos2SM + C * cosSigma * (-1.0 + 2.0 * cos2SM * cos2SM))); // (11)

      double delta = (lambda - lambdaOrig) / lambda;
      if (Math.abs(delta) < 1.0e-12) {
        break;
      }
    }

    double distance = (b * A * (sigma - deltaSigma));
    results[0] = distance;
    if (results.length > 1) {
      double initialBearing = Math.atan2(cosU2 * sinLambda, cosU1 * sinU2
          - sinU1 * cosU2 * cosLambda);
      initialBearing *= 180.0 / Math.PI;
      results[1] = initialBearing;
      if (results.length > 2) {
        double finalBearing = Math.atan2(cosU1 * sinLambda, -sinU1 * cosU2
            + cosU1 * sinU2 * cosLambda);
        finalBearing *= 180.0 / Math.PI;
        results[2] = finalBearing;
      }
    }
  }

  /*
   * Vincenty Direct Solution of Geodesics on the Ellipsoid (c) Chris Veness
   * 2005-2012
   * 
   * from: Vincenty direct formula - T Vincenty, "Direct and Inverse Solutions
   * of Geodesics on the Ellipsoid with application of nested equations", Survey
   * Review, vol XXII no 176, 1975 http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
   */

  /**
   * Calculates destination point and final bearing given given start point,
   * bearing & distance, using Vincenty inverse formula for ellipsoids
   * 
   * @param lat1
   *          start point latitude
   * @param lon1
   *          start point longitude
   * @param brng
   *          initial bearing in decimal degrees
   * @param dist
   *          distance along bearing in metres
   * @returns an array of the desination point coordinates and the final bearing
   */

  public static void computeDestinationAndBearing(double lat1, double lon1,
      double brng, double dist, double results[]) {
    double a = 6378137, b = 6356752.3142, f = 1 / 298.257223563; // WGS-84
                                                                 // ellipsiod
    double s = dist;
    double alpha1 = toRad(brng);
    double sinAlpha1 = Math.sin(alpha1);
    double cosAlpha1 = Math.cos(alpha1);

    double tanU1 = (1 - f) * Math.tan(toRad(lat1));
    double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1)), sinU1 = tanU1 * cosU1;
    double sigma1 = Math.atan2(tanU1, cosAlpha1);
    double sinAlpha = cosU1 * sinAlpha1;
    double cosSqAlpha = 1 - sinAlpha * sinAlpha;
    double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
    double A = 1 + uSq / 16384
        * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
    double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
    double sinSigma = 0, cosSigma = 0, deltaSigma = 0, cos2SigmaM = 0;
    double sigma = s / (b * A), sigmaP = 2 * Math.PI;

    while (Math.abs(sigma - sigmaP) > 1e-12) {
      cos2SigmaM = Math.cos(2 * sigma1 + sigma);
      sinSigma = Math.sin(sigma);
      cosSigma = Math.cos(sigma);
      deltaSigma = B
          * sinSigma
          * (cos2SigmaM + B
              / 4
              * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6
                  * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma)
                  * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
      sigmaP = sigma;
      sigma = s / (b * A) + deltaSigma;
    }

    double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
    double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
        (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
    double lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1
        * sinSigma * cosAlpha1);
    double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
    double L = lambda
        - (1 - C)
        * f
        * sinAlpha
        * (sigma + C * sinSigma
            * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    double lon2 = (toRad(lon1) + L + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalise
                                                                             // to
                                                                             // -180...+180

    double revAz = Math.atan2(sinAlpha, -tmp); // final bearing, if required

    results[0] = toDegrees(lat2);
    results[1] = toDegrees(lon2);
    results[2] = toDegrees(revAz);

  }

  private static double toRad(double angle) {
    return angle * Math.PI / 180;
  }

  private static double toDegrees(double radians) {
    return radians * 180 / Math.PI;
  }

}