Calculating coordinates given a bearing and a distance

Question

I am having problems implementing the function described here here.

This is my Java implementation:

private static double[] pointRadialDistance(double lat1, double lon1, 
        double radianBearing, double radialDistance) {
     double lat = Math.asin(Math.sin(lat1)*Math.cos(radialDistance)+Math.cos(lat1)
             *Math.sin(radialDistance)*Math.cos(radianBearing));
     double lon;
     if(Math.cos(lat) == 0) {  // Endpoint a pole
        lon=lon1;      
     }
     else {
        lon = ((lon1-Math.asin(Math.sin(radianBearing)*Math.sin(radialDistance)/Math.cos(lat))
                +Math.PI) % (2*Math.PI)) - Math.PI;
     }
    return (new double[]{lat, lon});
}

I convert the degree bearing to radians and convert the distance (km) into a radians distance before calling the function - so that's not the problem.

However, when I input coordinates such as: lat = 49.25705; lon = -123.140259; with a bearing of 225 (south-west) and a distance of 1km

I get this returned: lat: -1.0085434360125864 lon: -3.7595299668539504

Its obviously not correct, can anyone see what I am doing wrong?

Thanks

Answer 1

It seems like these are the issues in your code:

1. You need to convert lat1 and lon1 to radians before calling your function.
2. You may be scaling radialDistance incorrectly.
3. Testing a floating-point number for equality is dangerous. Two numbers that are equal after exact arithmetic might not be exactly equal after floating-point arithmetic. Thus abs(x-y) < threshold is safer than x == y for testing two floating-point numbers x and y for equality.
4. I think you want to convert lat and lon from radians to degrees.

Here is my implementation of your code in Python:

#!/usr/bin/env python

from math import asin,cos,pi,sin

rEarth = 6371.01 # Earth's average radius in km
epsilon = 0.000001 # threshold for floating-point equality


def deg2rad(angle):
    return angle*pi/180


def rad2deg(angle):
    return angle*180/pi


def pointRadialDistance(lat1, lon1, bearing, distance):
    """
    Return final coordinates (lat2,lon2) [in degrees] given initial coordinates
    (lat1,lon1) [in degrees] and a bearing [in degrees] and distance [in km]
    """
    rlat1 = deg2rad(lat1)
    rlon1 = deg2rad(lon1)
    rbearing = deg2rad(bearing)
    rdistance = distance / rEarth # normalize linear distance to radian angle

    rlat = asin( sin(rlat1) * cos(rdistance) + cos(rlat1) * sin(rdistance) * cos(rbearing) )

    if cos(rlat) == 0 or abs(cos(rlat)) < epsilon: # Endpoint a pole
        rlon=rlon1
    else:
        rlon = ( (rlon1 - asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi

    lat = rad2deg(rlat)
    lon = rad2deg(rlon)
    return (lat, lon)


def main():
    print "lat1 \t lon1 \t\t bear \t dist \t\t lat2 \t\t lon2"
    testcases = []
    testcases.append((0,0,0,1))
    testcases.append((0,0,90,1))
    testcases.append((0,0,0,100))
    testcases.append((0,0,90,100))
    testcases.append((49.25705,-123.140259,225,1))
    testcases.append((49.25705,-123.140259,225,100))
    testcases.append((49.25705,-123.140259,225,1000))
    for lat1, lon1, bear, dist in testcases:
        (lat,lon) = pointRadialDistance(lat1,lon1,bear,dist)
        print "%6.2f \t %6.2f \t %4.1f \t %6.1f \t %6.2f \t %6.2f" % (lat1,lon1,bear,dist,lat,lon)


if __name__ == "__main__":
    main()

Here is the output:

lat1     lon1        bear    dist        lat2        lon2
  0.00     0.00       0.0       1.0        0.01        0.00
  0.00     0.00      90.0       1.0        0.00       -0.01
  0.00     0.00       0.0     100.0        0.90        0.00
  0.00     0.00      90.0     100.0        0.00       -0.90
 49.26   -123.14     225.0      1.0       49.25      -123.13
 49.26   -123.14     225.0    100.0       48.62      -122.18
 49.26   -123.14     225.0   1000.0       42.55      -114.51

Answer 2

Fundamentally, it appears that your problem is that you are passing latitude, longitude and bearing as degrees rather than radians. Try ensuring that you are always passing radians to your function and see what you get back.

PS: see similar issues discussed here and here.

Answer 3

I think there is a problem in the Algorithm provided in message 5.

It works but for only for the latitude, for the longitude there is a problem because of the sign.

The data speaks for themselves :

49.26 -123.14 225.0 1.0 49.25 -123.13

If you start from -123.14° and go WEST you should have something FAR in the WEST. Here we go back on the EAST (-123.13) !

The formula should includes somewhere :

degreeBearing = ((360-degreeBearing)%360)

before radian convertion.

Answer 4

When I implemented this, my resulting latitudes were correct but the longitudes were wrong. For example starting point: 36.9460678N 9.434807E, Bearing 45.03334, Distance 15.0083313km The result was 37.0412865N 9.315302E That's further west than my starting point, rather than further east. In fact it's as if the bearing was 315.03334 degrees.

More web searching led me to: http://www.movable-type.co.uk/scripts/latlong.html The longitude code is show below (in C# with everything in radians)

        if ((Math.Cos(rLat2) == 0) || (Math.Abs(Math.Cos(rLat2)) < EPSILON))
        {
            rLon2 = rLon1;
        }
        else
        {
            rLon2 = rLon1 + Math.Atan2(Math.Sin(rBearing) * Math.Sin(rDistance) * Math.Cos(rLat1), Math.Cos(rDistance) - Math.Sin(rLat1) * Math.Sin(rLat2));
        }

This seems to work fine for me. Hope it's helpful.

Answer 5

Thanks for your python code I tried setting it up in my use case where I'm trying to find the lat lon of a point in between two others at a set distance from the first point so it's quite similare to your code appart that my bearing is dynamically calculated

startpoint(lat1) lon1/lat1 = 55.625541,-21.142463

end point (lat2) lon2/lat2 = 55.625792,-22.142248

my result should be a point in between these two at lon3/lat3 unfortunetly I get lon3/lat3 = 0.0267695450609,0.0223553243666

I thought this might be a difference in lat lon but no when I add or sub it it's not good

any advice would be really great Thanks

here's my implementation

distance = 0.001 epsilon = 0.000001

calculating bearing dynamically

y = math.sin(distance) * math.cos(lat2);
x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(distance);
bearing = math.atan2(y, x)

calculating lat3 lon3 dynamically

rlat1 = (lat1 * 180) / math.pi
rlon1 = (lon1 * 180) / math.pi
rbearing = (bearing * 180) / math.pi
rdistance = distance / R # normalize linear distance to radian angle

rlat = math.asin( math.sin(rlat1) * math.cos(rdistance) + math.cos(rlat1) * math.sin(rdistance) * math.cos(rbearing) )
if math.cos(rlat) == 0 or abs(math.cos(rlat)) < epsilon: # Endpoint a pole
      rlon=rlon1
else:
    rlon = ( (rlon1 + math.asin( math.sin(rbearing)* math.sin(rdistance) / math.cos(rlat) ) + math.pi ) % (2*math.pi) ) - math.pi

lat3 = (rlat * math.pi)/ 180
lon3 = (rlon * math.pi)/ 180

Answer 6

Everything is working as intended, but the problem is that your maths assumes the Earth to be a sphere when in reality it approximates an ellipsoid.

A quick trawl of your favoured search engine for 'Vincenty Formula' will hopefully prove useful.