CODEWITHSUNDEEP
pythonfunctionlatitude-longitudeplotly-dash
Sterling Butters
Sterling Butters

Reputation: 1175

Python: Calculate bearing between two lat/long

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.

Upvotes: 20

Views: 40039

Answers (5)

Al Zonke
Al Zonke

Reputation: 552

geopy uses geographiclib as geodesic engine. It has all needed functionality. Refrain from using custom trigonometric implementation — is way to inaccurate.

from geographiclib.geodesic import Geodesic

def get_bearing(point_a: Point, point_b: Point) -> float:
    result = Geodesic.WGS84.Inverse(point_a.latitude, point_a.longitude, point_b.latitude, point_b.longitude,)
    bearing = result['azi1'] % 360
    return bearing

Upvotes: 0

J. Taylor
J. Taylor

Reputation: 4855

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.

Upvotes: 49

aliff danial
aliff danial

Reputation: 119

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

Upvotes: 11

Omkar Mestry
Omkar Mestry

Reputation: 41

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/

Upvotes: 2

Sterling Butters
Sterling Butters

Reputation: 1175

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

Upvotes: 11

Related Questions