CODEWITHSUNDEEP
pythongpscoordinatesbearing
242Eld
242Eld

Reputation: 225

calculating a gps coordinate given a point, bearing and distance

I have a problem which draws my back in some project for some time now.

I'm basically looking to trap a polygon using x, y points drawn by some script I've written. lat, lon are the center GPS cords of the polygon and I'm looking for its surrounding polygon.

here is a part of my code in python:

def getcords(lat, lon, dr, bearing):
    lat2=asin(sin(lat)*cos(dr)+cos(lat)*sin(dr)*cos(bearing))
    lon2=lon+atan2(sin(bearing)*sin(dr)*cos(lat),cos(dr)-sin(lat)*sin(lat2))
    return [lat2,lon2]

my input goes like this:

However for the input:

getcorsds1(42.189275, -76.85823, 0.5/3958.82, 30)

I get output: [-1.3485899508698462, -76.8576637627568], however [42.2516666666667, -76.8097222222222] is the right answer.

as for the angular distance, I calculate it simply by dividing the distance in miles by the earth's radius(=3958.82).

anybody?

Upvotes: 5

Views: 15643

Answers (3)

tacosan
tacosan

Reputation: 194

With geopy v2.0.0 (+ kilometers instead miles)

from geopy import Point                                                                                                                                                                       
from geopy.distance import geodesic                                                                                                                                                           
                                                                                                                                                                                              
distKm = 1                                                                                                                                                                                    
lat1 = 35.68096477080332                                                                                                                                                                      
lon1 = 139.76720809936523                                                                                                                                                                     
                                                                                                                                                                                              
print('center', lat1, lon1)                                                                                                                                                                   
print('north', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 0).format_decimal())                                                                                                
print('east', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 90).format_decimal())                                                                                                
print('south', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 180).format_decimal())                                                                                              
print('west', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 270).format_decimal())

result is

center 35.6809647708 139.767208099
north 35.6899775841, 139.767208099
east 35.680964264, 139.778254714
south 35.6719519439, 139.767208099
west 35.680964264, 139.756161485

Upvotes: 5

John Machin
John Machin

Reputation: 83032

The sin and cos functions expect their arguments in radians, not in degrees. The asin and atan2 functions produce a result in radians, not in degrees. In general, one needs to convert input angles (lat1, lon1 and bearing) from degrees to radians using math.radians() and convert output angles (lat2 and lon2) from radians to degrees using math.degrees().

Note that your code has two other problems:

(1) It doesn't allow for travel across the 180-degrees meridian of longitude; you need to constrain your answer such that -180 <= longitude_degrees <= +180.

(2) If you are going to use this function extensively, you might like to remove the redundant calculations: sin(lat1), cos(dr), cos(lat1), and sin(dr) are each calculated twice.

Upvotes: 3

eumiro
eumiro

Reputation: 213125

Why don't you use nice libraries?

from geopy import Point
from geopy.distance import distance, VincentyDistance

# given: lat1, lon1, bearing, distMiles
lat2, lon2 = VincentyDistance(miles=distMiles).destination(Point(lat1, lon1), bearing)

For lat1, lon1, distMiles, bearing = 42.189275,-76.85823, 0.5, 30 it returns 42.1955489, -76.853359.

Upvotes: 4

Related Questions