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

Solution 1

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

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 = math.radians(52.20472) #Current lat point converted to radians
lon1 = math.radians(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.degrees(lat2)
lon2 = math.degrees(lon2)

print(lat2)
print(lon2)

Solution 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

Solution 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

Solution 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;
 }

Solution 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:

View more solutions

Comments

[Peter]

Not working. From North to South, result is correct but it's wrong in "East-West" direction.

[user2360915]

Look good, but i think the requestor would like to have something in python. Wrong?

[jpoehnelt]

this library has some distance issues waiting to be resolved: github.com/geopy/geopy/pull/144

[Dev Null]

may be better named get_gps_coord or similar. You're not getting the distance, you supply that to the func. But thanks for this, it's exactly what I was looking for. Many searches return calculating distance between coords (false positives). Thanks!

[not2qubit]

Please post functional code, including what it need to run. I.e. this seem to be dependent on Microsoft.Maps. Where to find/ how to install this?

[Brad Mathews]

You would only use Bing (Microsoft) Maps if your program uses Bing maps. Just take the Math.degrees(lat) and Math.degrees(lon) values and do with them whatever you need to for your application.

[boden]

same result for me as well

[Tommie C.]

Could you explain a couple of the variables? $range and $magvar could use a bit more exposition for the novice readers like (me:)

[cloudscomputes]

The geojson's map is very helpful for me to target a location in a map

[jkamizato]

Awesome! Thanks for your contribution!

[Vsevolod Krasnov]

Thank you, implemented that snippet in Kotlin.

[not2qubit]

Without stating the method you're using for calculation, the answer is basically useless.

[not2qubit]

Please see my answer and link to the formula that it uses and to what accuracy we can expect.

[hlongmore]

@not2qubit Whether @plinio-bueno-andrade-silva was aware or not, geopy.distance.distance currently uses geodesic. geopy And to be more specific, the ellipsoidal model used by default is the WGS-84 ellipsoid, "which is the most globally accurate."

[JVE999]

6,378.14 km seems to be the maximum radius of Earth. The average is about 6,371.0 km, which may allow for more accurate calculations.

[JVE999]

I noticed that if the original latitude is 0, the original longitude is -179, the bearing is 270 degrees (1.5pi radians), and the distance is 1500km, the resulting longitude is -192.4, which doesn't exist on a map.

[Mark Wardell]

Thank you implemented a snippet in C# gist.github.com/BicycleMark/3e1a2152febaa2935e4c8cfcea7e061b

[Michael Behrens]

I validated the code output using: fcc.gov/media/radio/find-terminal-coordinates

[Douglas Schmidt]

Thank you @kodisha, your fiddle help me a lot!

[Robert Moskal]

Thanks for saving me a little time.

[danjampro]

Note that the API has changed since v2.0.0. Instead use geopy.distance.geodesic: stackoverflow.com/a/62866744/4717384

[afp]

I think the last part of the longitude computation might be wrong, since variable lat is already updated before calculating lon, i.e. the term sin(lat) * sin(lat) is not actually using both the old and the new latitudes, respectively.

[afp]

Same as my comment in the previous answer, I think the last part of the longitude computation might be wrong, since variable lat is already updated before calculating lon, i.e. the term Math.sin(lat) * Math.sin(lat) is not actually using both the old and the new latitudes, respectively.

[afp]

Probably this is the most correct answer, since it correctly uses the old and new latitudes, respectively, when calculating the last term of the lon2 expression, i.e. Math.sin(lat1)*Math.sin(lat2). Hence the slightly different result.

[afp]

Most answers don't work properly because of a small mistake in the last term; they don't use the old & new latitudes, but the already updated latitude twice.