20

I am attempting to calculate the bearing between two lat/long.

I don't have a question regarding the function/formula per se,

provided:

def get_bearing(lat1, long1, lat2, long2):
    dLon = (long2 - long1)

    y = math.sin(dLon) * math.cos(lat2)
    x = math.cos(lat1) * math.sin(lat2) - math.sin(lat1) * math.cos(lat2) * math.cos(dLon)

    brng = math.atan2(y, x)

    brng = np.rad2deg(brng)

    return brng

the problem is that the result isn't what is expected.

The intended usage of the function returns the bearing between two lat/long pairs in a (very long) list i.e.

    lat1 = path[int(len(path) * location / 1000)][0]
    lat2 = path[int(len(path) * location / 1000) + 1][0]
    lng1 = path[int(len(path) * location / 1000)][1]
    lng2 = path[int(len(path) * location / 1000) + 1][1]

The bearing result then alters the view orientation of the plot where bearing can assume a value in the range [-180, 180]. Ideally, the result would appear such that the line formed between lat1, lng1 and lat2, lng2 is perfectly "vertical" in the plot (lat/lon annotations are switched in plot), see below

enter image description here

enter image description here

I am hoping that someone might be able to deduce the problem from the bearing returned from the function and what the expected bearing should be. A few instances below:

Current Location: 30.07134 -97.23076
Next in path: 30.0709 -97.22907
Calculated Bearing: 88.39967863143139
Expected Bearing: ~-70.67

Current Location: 29.91581 -96.85068
Next in path: 29.91556 -96.85021
Calculated Bearing: 118.9170342272798
Expected Bearing: ~122.67

Current Location: 29.69419 -96.53487
Next in path: 29.69432 -96.53466
Calculated Bearing 141.0271357781952
Expected Bearing: ~56

Current Location: 29.77357 -96.07924
Next in path: 29.77349 -96.07876
Calculated Bearing 165.24612555483893
Expected Bearing: ~104

Happy to provide additional information, thanks in advance for any/all help.

4 Answers 4

48

Have you considered using pyproj to do the calculations instead of rolling your own?:

import pyproj
geodesic = pyproj.Geod(ellps='WGS84')
fwd_azimuth,back_azimuth,distance = geodesic.inv(long1, lat1, long2, lat2)

In this example fwd_azimuth is the bearing you are after and back_azimuth is inverse bearing (going the opposite direction).

I used WGS84 here, so you would need to replace with correct coordinate system, and need to rewrite ensure lat/long are correct type of coordinates for geodesic.inv(). But using a well-tested, existing geo-spatial lib will likely save you a lot of hair pulling.

5
  • 1
    I'm getting a nan returned for fwd_azimuth and that is with float(lat), float(lon)..., etc - I would definitely like to use a geo-spatial lib but I didn't feel like dealing with idiosyncratic issues, any ideas on the nan? Commented Feb 25, 2019 at 20:47
  • What coordinate system are you using for your input lat/long pairs? Take a look at the API docs for pyproj to see how they represent lat/long. Like I said above, you'd have to modify your code to pass in lat/long in this format.
    – J. Taylor
    Commented Feb 25, 2019 at 20:51
  • tbh, idk - I am using Plot.ly's ScatterMapBox object and I can't find anything in their reference: plot.ly/python/reference/#scattermapbox. I'm going to guess it is the same as GoogleMaps given that their API works well with ScatterMapBox. Google Maps and Microsoft Virtual Earth use a Mercator projection based on the World Geodetic System (WGS) 1984 geographic coordinate system (datum). Commented Feb 25, 2019 at 20:56
  • 3
    I found a lib (geographiclib) that I like and that works! I will upvote your answer for leading me in that direction Commented Feb 25, 2019 at 21:52
  • 3
    You have lats and longs the wrong way around, see pyproj4.github.io/pyproj/stable/api/geod.html in there for example: >>> from pyproj import Geod >>> g = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid. >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.) >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.) >>> # compute forward and back azimuths, plus distance >>> # between Boston and Portland. >>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat) >>> "%7.3f %6.3f %12.3f" % (az12,az21,dist) '-66.531 75.654 4164192.708' Commented Jul 9, 2020 at 14:23
11

Ended up changing the function:

from geographiclib.geodesic import Geodesic
...
def get_bearing(lat1, lat2, long1, long2):
    brng = Geodesic.WGS84.Inverse(lat1, long1, lat2, long2)['azi1']
    return brng
1
10

The code for finding the bearing is fine. But you just need to add math.radians when finding X and Y.

import numpy
import math

def get_bearing(lat1, long1, lat2, long2):
    dLon = (long2 - long1)
    x = math.cos(math.radians(lat2)) * math.sin(math.radians(dLon))
    y = math.cos(math.radians(lat1)) * math.sin(math.radians(lat2)) - math.sin(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.cos(math.radians(dLon))
    brng = numpy.arctan2(x,y)
    brng = numpy.degrees(brng)

    return brng
2
  • The code snippet is missing imports - where are arctan2 and degrees imported from, should that be numpy? Commented Jul 13, 2021 at 13:36
  • Yeah it is imported from numpy. Sorry I didnt put it there. @Wolkenarchitekt Commented Jul 16, 2021 at 8:25
2

You can try this new package called Turfpy which is a porting of turf.js in python.

from turfpy import measurement
from geojson import Point, Feature
start = Feature(geometry=Point((-75.343, 39.984)))
end = Feature(geometry=Point((-75.534, 39.123)))
measurement.bearing(start,end)

https://pypi.org/project/turfpy/

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.