Get lat/long given current point, distance and bearing

Question

Given an existing point in lat/long, distance in (in KM) and bearing (in degrees converted to radians), I would like to calculate the new lat/long. This site crops up over and over again, but I just can't get the formula to work for me.

The formulas as taken the above link are:

lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))

lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))

The above formula is for MSExcel where-

asin          = arc sin()   
d             = distance (in any unit)   
R             = Radius of the earth (in the same unit as above)  
and hence d/r = is the angular distance (in radians)  
atan2(a,b)    = arc tan(b/a)  
θ is the bearing (in radians, clockwise from north);  

Here's the code I've got in Python.

import math

R = 6378.1 #Radius of the Earth
brng = 1.57 #Bearing is 90 degrees converted to radians.
d = 15 #Distance in km

#lat2  52.20444 - the lat result I'm hoping for
#lon2  0.36056 - the long result I'm hoping for.

lat1 = 52.20472 * (math.pi * 180) #Current lat point converted to radians
lon1 = 0.14056 * (math.pi * 180) #Current long point converted to radians

lat2 = math.asin( math.sin(lat1)*math.cos(d/R) +
             math.cos(lat1)*math.sin(d/R)*math.cos(brng))

lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1),
                     math.cos(d/R)-math.sin(lat1)*math.sin(lat2))

print(lat2)
print(lon2)

I get

lat2 = 0.472492248844 
lon2 = 79.4821662373

Answer 1

Needed to convert answers from radians back to degrees. Working code below:

from math import asin, atan2, cos, degrees, radians, sin

def get_point_at_distance(lat1, lon1, d, bearing, R=6371):
    """
    lat: initial latitude, in degrees
    lon: initial longitude, in degrees
    d: target distance from initial
    bearing: (true) heading in degrees
    R: optional radius of sphere, defaults to mean radius of earth

    Returns new lat/lon coordinate {d}km from initial, in degrees
    """
    lat1 = radians(lat1)
    lon1 = radians(lon1)
    a = radians(bearing)
    lat2 = asin(sin(lat1) * cos(d/R) + cos(lat1) * sin(d/R) * cos(a))
    lon2 = lon1 + atan2(
        sin(a) * sin(d/R) * cos(lat1),
        cos(d/R) - sin(lat1) * sin(lat2)
    )
    return (degrees(lat2), degrees(lon2),)

lat = 52.20472
lon = 0.14056
distance = 15
bearing = 90
lat2, lon2 = get_point_at_distance(lat, lon, distance, bearing)

# lat2  52.20444 - the lat result I'm hoping for
# lon2  0.36056 - the long result I'm hoping for.

print(lat2, lon2)
# prints "52.20451523755824 0.36067845713550956"

Answer 2

The geopy library supports this:

import geopy
from geopy.distance import VincentyDistance

# given: lat1, lon1, b = bearing in degrees, d = distance in kilometers

origin = geopy.Point(lat1, lon1)
destination = VincentyDistance(kilometers=d).destination(origin, b)

lat2, lon2 = destination.latitude, destination.longitude

Found via https://stackoverflow.com/a/4531227/37610

Answer 3

This question is known as the direct problem in the study of geodesy.

This is indeed a very popular question and one that is a constant cause of confusion. The reason is that most people are looking for a simple and straight-forward answer. But there is none, because most people asking this question are not supplying enough information, simply because they are not aware that:

1. Earth is not a perfect sphere, since it is flattened/compressed by it poles
2. Because of (1) earth does not have a constant Radius, R. See here.
3. Earth is not perfectly smooth (variations in altitude) etc.
4. Due to tectonic plate movement, a geographic point's lat/lon position may change by several millimeters (at least), every year.

Therefore there are many different assumptions used in the various geometric models that apply differently, depending on your needed accuracy. So to answer the question you need to consider to what accuracy you would like to have your result.

Some examples:

• I'm just looking for an approximate location to the nearest few kilometers for small ( < 100 km) distances of in latitudes between 0-70 deg N|S. (Earth is ~flat model.)
• I want an answer that is good anywhere on the globe, but only accurate to about a few meters
• I want a super accurate positioning that is valid down to atomic scales of nanometers [nm].
• I want answers that is very fast and easy to calculate and not computationally intensive.

So you can have many choices in which algorithm to use. In addition each programming language has it's own implementation or "package" multiplied by number of models and the model developers specific needs. For all practical purposes here, it pays off to ignore any other language apart javascript, since it very closely resemble pseudo-code by its nature. Thus it can be easily converted to any other language, with minimal changes.

Then the main models are:

• Euclidian/Flat earth model: good for very short distances under ~10 km
• Spherical model: good for large longitudinal distances, but with small latitudinal difference. Popular model:
  • Haversine: meter accuracy on [km] scales, very simple code.
• Ellipsoidal models: Most accurate at any lat/lon and distance, but is still a numerical approximation that depend on what accuracy you need. Some popular models are:
  • Lambert: ~10 meter precision over 1000's of km.
  • Paul D.Thomas: Andoyer-Lambert approximation
  • Vincenty: millimeter precision and computational efficiency
  • Kerney: nanometer precision

References:

• https://en.wikipedia.org/wiki/Reference_ellipsoid
• https://en.wikipedia.org/wiki/Haversine_formula
• https://en.wikipedia.org/wiki/Earth_ellipsoid
• https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid
• https://en.wikipedia.org/wiki/Vincenty%27s_formulae
• https://geographiclib.sourceforge.io/scripts/geod-calc.html

Answer 4

May be a bit late for answering, but after testing the other answers, it appears they don't work correctly. Here is a PHP code we use for our system. Working in all directions.

PHP code:

lat1 = latitude of start point in degrees

long1 = longitude of start point in degrees

d = distance in KM

angle = bearing in degrees

function get_gps_distance($lat1,$long1,$d,$angle)
{
    # Earth Radious in KM
    $R = 6378.14;

    # Degree to Radian
    $latitude1 = $lat1 * (M_PI/180);
    $longitude1 = $long1 * (M_PI/180);
    $brng = $angle * (M_PI/180);

    $latitude2 = asin(sin($latitude1)*cos($d/$R) + cos($latitude1)*sin($d/$R)*cos($brng));
    $longitude2 = $longitude1 + atan2(sin($brng)*sin($d/$R)*cos($latitude1),cos($d/$R)-sin($latitude1)*sin($latitude2));

    # back to degrees
    $latitude2 = $latitude2 * (180/M_PI);
    $longitude2 = $longitude2 * (180/M_PI);

    # 6 decimal for Leaflet and other system compatibility
   $lat2 = round ($latitude2,6);
   $long2 = round ($longitude2,6);

   // Push in array and get back
   $tab[0] = $lat2;
   $tab[1] = $long2;
   return $tab;
 }

Answer 5

I ported answer by Brad to vanilla JS answer, with no Bing maps dependency

https://jsfiddle.net/kodisha/8a3hcjtd/

    // ----------------------------------------
    // Calculate new Lat/Lng from original points
    // on a distance and bearing (angle)
    // ----------------------------------------
    let llFromDistance = function(latitude, longitude, distance, bearing) {
      // taken from: https://stackoverflow.com/a/46410871/13549 
      // distance in KM, bearing in degrees
    
      const R = 6378.1; // Radius of the Earth
      const brng = bearing * Math.PI / 180; // Convert bearing to radian
      let lat = latitude * Math.PI / 180; // Current coords to radians
      let lon = longitude * Math.PI / 180;
    
      // Do the math magic
      lat = Math.asin(Math.sin(lat) * Math.cos(distance / R) + Math.cos(lat) * Math.sin(distance / R) * Math.cos(brng));
      lon += Math.atan2(Math.sin(brng) * Math.sin(distance / R) * Math.cos(lat), Math.cos(distance / R) - Math.sin(lat) * Math.sin(lat));
    
      // Coords back to degrees and return
      return [(lat * 180 / Math.PI), (lon * 180 / Math.PI)];
    
    }
    
    let pointsOnMapCircle = function(latitude, longitude, distance, numPoints) {
      const points = [];
      for (let i = 0; i <= numPoints - 1; i++) {
        const bearing = Math.round((360 / numPoints) * i);
        console.log(bearing, i);
        const newPoints = llFromDistance(latitude, longitude, distance, bearing);
        points.push(newPoints);
      }
      return points;
    }
    
    const points = pointsOnMapCircle(41.890242042122836, 12.492358982563019, 0.2, 8);
    let geoJSON = {
      "type": "FeatureCollection",
      "features": []
    };
    points.forEach((p) => {
      geoJSON.features.push({
        "type": "Feature",
        "properties": {},
        "geometry": {
          "type": "Point",
          "coordinates": [
            p[1],
            p[0]
          ]
        }
      });
    });
    
    document.getElementById('res').innerHTML = JSON.stringify(geoJSON, true, 2);

In addition, I added geoJSON export, so you can simply paste resulting geoJSON to: http://geojson.io/#map=17/41.89017/12.49171 to see the results instantly.

Result:

Answer 6

Quick way using geopy

from geopy import distance
#distance.distance(unit=15).destination((lat,lon),bering) 
#Exemples
distance.distance(nautical=15).destination((-24,-42),90) 
distance.distance(miles=15).destination((-24,-42),90)
distance.distance(kilometers=15).destination((-24,-42),90) 

Answer 7

lon1 and lat1 in degrees

brng = bearing in radians

d = distance in km

R = radius of the Earth in km

lat2 = math.degrees((d/R) * math.cos(brng)) + lat1
long2 = math.degrees((d/(R*math.sin(math.radians(lat2)))) * math.sin(brng)) + long1

I implemented your algorithm and mine in PHP and benchmarked it. This version ran in about 50% of the time. The results generated were identical, so it seems to be mathematically equivalent.

I didn't test the python code above so there might be syntax errors.

Answer 8

I ported the Python to Javascript. This returns a Bing Maps Location object, you can change to whatever you like.

getLocationXDistanceFromLocation: function(latitude, longitude, distance, bearing) {
    // distance in KM, bearing in degrees

    var R = 6378.1,                         // Radius of the Earth
        brng = Math.radians(bearing)       // Convert bearing to radian
        lat = Math.radians(latitude),       // Current coords to radians
        lon = Math.radians(longitude);

    // Do the math magic
    lat = Math.asin(Math.sin(lat) * Math.cos(distance / R) + Math.cos(lat) * Math.sin(distance / R) * Math.cos(brng));
    lon += Math.atan2(Math.sin(brng) * Math.sin(distance / R) * Math.cos(lat), Math.cos(distance/R)-Math.sin(lat)*Math.sin(lat));

    // Coords back to degrees and return
    return new Microsoft.Maps.Location(Math.degrees(lat), Math.degrees(lon));

},

Answer 9

Thanks to @kodisha, here is a Swift version, but with improved and more precise calculation for Earth radius:

extension CLLocationCoordinate2D {
  
  func earthRadius() -> CLLocationDistance {
    let earthRadiusInMetersAtSeaLevel = 6378137.0
    let earthRadiusInMetersAtPole = 6356752.314
    
    let r1 = earthRadiusInMetersAtSeaLevel
    let r2 = earthRadiusInMetersAtPole
    let beta = latitude
    
    let earthRadiuseAtGivenLatitude = (
      ( pow(pow(r1, 2) * cos(beta), 2) + pow(pow(r2, 2) * sin(beta), 2) ) /
        ( pow(r1 * cos(beta), 2) + pow(r2 * sin(beta), 2) )
    )
    .squareRoot()
    
    return earthRadiuseAtGivenLatitude
  }
  
  func locationByAdding(
    distance: CLLocationDistance,
    bearing: CLLocationDegrees
  ) -> CLLocationCoordinate2D {
    let latitude = self.latitude
    let longitude = self.longitude
    
    let earthRadiusInMeters = self.earthRadius()
    let brng = bearing.degreesToRadians
    var lat = latitude.degreesToRadians
    var lon = longitude.degreesToRadians
    
    lat = asin(
      sin(lat) * cos(distance / earthRadiusInMeters) +
        cos(lat) * sin(distance / earthRadiusInMeters) * cos(brng)
    )
    lon += atan2(
      sin(brng) * sin(distance / earthRadiusInMeters) * cos(lat),
      cos(distance / earthRadiusInMeters) - sin(lat) * sin(lat)
    )
    
    let newCoordinate = CLLocationCoordinate2D(
      latitude: lat.radiansToDegrees,
      longitude: lon.radiansToDegrees
    )
    
    return newCoordinate
  }
}

extension FloatingPoint {
  var degreesToRadians: Self { self * .pi / 180 }
  var radiansToDegrees: Self { self * 180 / .pi }
}

Answer 10

Also late but for those who might find this, you will get more accurate results using the geographiclib library. Check out the geodesic problem descriptions and the JavaScript examples for an easy introduction to how to use to answer the subject question as well as many others. Implementations in a variety of languages including Python. Far better than coding your own if you care about accuracy; better than VincentyDistance in the earlier "use a library" recommendation. As the documentation says: "The emphasis is on returning accurate results with errors close to round-off (about 5–15 nanometers)."

Answer 11

Just interchange the values in the atan2(y,x) function. Not atan2(x,y)!

Answer 12

I ported the answer from @David M to java if anyone wanted this... I do get a slight different result of 52.20462299620793, 0.360433887489931

    double R = 6378.1;  //Radius of the Earth
    double brng = 1.57;  //Bearing is 90 degrees converted to radians.
    double d = 15;  //Distance in km

    double lat2 = 52.20444; // - the lat result I'm hoping for
    double lon2 = 0.36056; // - the long result I'm hoping for.

    double lat1 = Math.toRadians(52.20472); //Current lat point converted to radians
    double lon1 = Math.toRadians(0.14056); //Current long point converted to radians

    lat2 = Math.asin( Math.sin(lat1)*Math.cos(d/R) +
            Math.cos(lat1)*Math.sin(d/R)*Math.cos(brng));

    lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(d/R)*Math.cos(lat1),
            Math.cos(d/R)-Math.sin(lat1)*Math.sin(lat2));

    lat2 = Math.toDegrees(lat2);
    lon2 = Math.toDegrees(lon2);

    System.out.println(lat2 + ", " + lon2);

Answer 13

Here is a PHP version based on Ed Williams Aviation Formulary. Modulus is handled a little different in PHP. This works for me.

function get_new_waypoint ( $lat, $lon, $radial, $magvar, $range )
{

   // $range in nm.
   // $radial is heading to or bearing from    
   // $magvar for local area.

   $range = $range * pi() /(180*60);
   $radial = $radial - $magvar ;

   if ( $radial < 1 )
     {
        $radial = 360 + $radial - $magvar; 
     }
   $radial =  deg2rad($radial);
   $tmp_lat = deg2rad($lat);
   $tmp_lon = deg2rad($lon);
   $new_lat = asin(sin($tmp_lat)* cos($range) + cos($tmp_lat) * sin($range) * cos($radial));
   $new_lat = rad2deg($new_lat);
   $new_lon = $tmp_lon - asin(sin($radial) * sin($range)/cos($new_lat))+ pi() % 2 * pi() -  pi();
   $new_lon = rad2deg($new_lon);

   return $new_lat." ".$new_lon;

}

Answer 14

For whoever is interested in a Java solution here is my code: I noticed that the initial solution needs some tweaks in order to return a proper longitude value, especially when the point is at one of the poles. Also a round operation is sometimes required as the results on 0 latitude / longitude seem to slightly shift away from 0. For small distances, rounding will help in this regard.

private static final double EARTH_RADIUS = 6371; // average earth radius

    /**
     * Returns the coordinates of the point situated at the distance specified, in
     * the direction specified. Note that the value is an approximation, not an
     * exact result.
     *
     * @param startPointLatitude
     * @param startPointLongitude
     * @param distanceInKm
     * @param bearing:            0 means moving north, 90 moving east, 180 moving
     *                            south, 270 moving west. Max value 360 min value 0;
     * @return new point location
     */
    public static LocationDTO getPointAt(double startPointLatitude, double startPointLongitude, double distanceInKm,
            double bearing) {
        if (Math.abs(startPointLatitude) > 90) {
            throw new BadRequestException(ExceptionMessages.INVALID_LATITUDE);
        } else if (Math.abs(startPointLatitude) == 90) {
            startPointLatitude = 89.99999 * Math.signum(startPointLatitude); // we have to do this conversion else the formula doesnt return the correct longitude value 
        }
        if (Math.abs(startPointLongitude) > 180) {
            throw new BadRequestException(ExceptionMessages.INVALID_LONGITUDE);
        }
        double angularDistance = distanceInKm / EARTH_RADIUS;
        bearing = deg2rad(bearing);
        startPointLatitude = deg2rad(startPointLatitude);
        startPointLongitude = deg2rad(startPointLongitude);
        double latitude = Math.asin(Math.sin(startPointLatitude) * Math.cos(angularDistance)
                + Math.cos(startPointLatitude) * Math.sin(angularDistance) * Math.cos(bearing));
        double longitude = startPointLongitude
                + Math.atan2(Math.sin(bearing) * Math.sin(angularDistance) * Math.cos(startPointLatitude),
                        Math.cos(angularDistance) - Math.sin(startPointLatitude) * Math.sin(latitude));
        longitude = (rad2deg(longitude) + 540) % 360 - 180; // normalize longitude to be in -180 +180 interval 
        LocationDTO result = new LocationDTO();
        result.setLatitude(roundValue(rad2deg(latitude)));
        result.setLongitude(roundValue(longitude));
        return result;
    }

    private static double roundValue(double value) {
        DecimalFormat df = new DecimalFormat("#.#####");
        df.setRoundingMode(RoundingMode.CEILING);
        return Double.valueOf(df.format(value));
    }


    // This function converts decimal degrees to radians
    private static double deg2rad(double deg) {
        return (deg * Math.PI / 180.0);
    }

    // This function converts radians to decimal degrees
    private static double rad2deg(double rad) {
        return (rad * 180.0 / Math.PI);
    }

Answer 15

Very late to the party, but here is answer in R for anyone interested. Only change I've made is that I've set the radius to metres, so d needs to be set to meters too.

get_point_at_distance <- function(lon, lat, d, bearing, R = 6378137) {
  # lat: initial latitude, in degrees
  # lon: initial longitude, in degrees
  # d: target distance from initial point (in m)
  # bearing: (true) heading in degrees
  # R: mean radius of earth (in m)
  # Returns new lat/lon coordinate {d} m from initial, in degrees
  ## convert to radians
  lat1 <- lat * (pi/180)
  lon1 <- lon * (pi/180)
  a <- bearing * (pi/180)
  ## new position
  lat2 <- asin(sin(lat1) * cos(d/R) + cos(lat1) * sin(d/R) * cos(a))
  lon2 <- lon1 + atan2(
    sin(a) * sin(d/R) * cos(lat1),
    cos(d/R) - sin(lat1) * sin(lat2)
  )
  ## convert back to degrees
  lat2 <- lat2 * (180/pi)
  lon2 <- lon2 * (180/pi)
  ## return
  return(c(lon2, lat2))
}

lat = 52.20472
lon = 0.14056
distance = 15000
bearing = 90
get_point_at_distance(lon = lon, lat = lat, d = distance, bearing = bearing)
# [1]  0.3604322 52.2045157