CODEWITHSUNDEEP
algorithmgpscoordinatesbearing
user106996
user106996

Reputation: 191

Calculating coordinates given a bearing and a distance

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

Upvotes: 6

Views: 36981

Answers (7)

Vishal Gaddam
Vishal Gaddam

Reputation: 1

# -*- coding: utf-8 -*-
from math import asin, cos, pi, sin
import pandas as pd`enter code here`

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):
    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():
    source_path = "Source_File.xlsx"
    destination_path = "Destination_File.xlsx"
    
    df = pd.read_excel(source_path)
    for index, row in df.iterrows():
        latitude = row['Latitude']
        longitude = row['Longitude']
        bearing = 60
        distance = 0.060
        (lat,lon) = pointRadialDistance(latitude, longitude, bearing, distance)
        df.loc[index, 'New Latitude'] = lat
        df.loc[index, 'New Longitude'] = lon
        
    df.to_excel(destination_path, index=False)


if __name__ == "__main__":
    main()

Upvotes: 0

las3rjock
las3rjock

Reputation: 8724

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

Upvotes: 19

St&#233;phane
St&#233;phane

Reputation: 41

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.

Upvotes: 4

Ximines
Ximines

Reputation: 9

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.

Upvotes: 0

Tibor
Tibor

Reputation: 41

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

Upvotes: 0

user164104
user164104

Reputation:

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.

Upvotes: 1

Bob Cross
Bob Cross

Reputation: 22292

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.

Upvotes: 3

Related Questions