Haversine Formula in Python (Bearing and Distance between two GPS points)

Question

Problem

I would like to know how to get the distance and bearing between 2 GPS points. I have researched on the haversine formula. Someone told me that I could also find the bearing using the same data.

Edit

Everything is working fine but the bearing doesn't quite work right yet. The bearing outputs negative but should be between 0 - 360 degrees. The set data should make the horizontal bearing 96.02166666666666 and is:

Start point: 53.32055555555556 , -1.7297222222222221   
Bearing:  96.02166666666666  
Distance: 2 km  
Destination point: 53.31861111111111, -1.6997222222222223  
Final bearing: 96.04555555555555

Here is my new code:

from math import *

Aaltitude = 2000
Oppsite  = 20000

lat1 = 53.32055555555556
lat2 = 53.31861111111111
lon1 = -1.7297222222222221
lon2 = -1.6997222222222223

lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * atan2(sqrt(a), sqrt(1-a))
Base = 6371 * c


Bearing =atan2(cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1), sin(lon2-lon1)*cos(lat2)) 

Bearing = degrees(Bearing)
print ""
print ""
print "--------------------"
print "Horizontal Distance:"
print Base
print "--------------------"
print "Bearing:"
print Bearing
print "--------------------"


Base2 = Base * 1000
distance = Base * 2 + Oppsite * 2 / 2
Caltitude = Oppsite - Aaltitude

a = Oppsite/Base
b = atan(a)
c = degrees(b)

distance = distance / 1000

print "The degree of vertical angle is:"
print c
print "--------------------"
print "The distance between the Balloon GPS and the Antenna GPS is:"
print distance
print "--------------------"

Solution 1

Here's a Python version:

from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance in kilometers between two points 
    on the earth (specified in decimal degrees)
    """
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a)) 
    r = 6371 # Radius of earth in kilometers. Use 3956 for miles. Determines return value units.
    return c * r

Solution 2

Most of these answers are "rounding" the radius of the earth. If you check these against other distance calculators (such as geopy), these functions will be off.

This works well:

from math import radians, cos, sin, asin, sqrt

def haversine(lat1, lon1, lat2, lon2):

      R = 3959.87433 # this is in miles.  For Earth radius in kilometers use 6372.8 km

      dLat = radians(lat2 - lat1)
      dLon = radians(lon2 - lon1)
      lat1 = radians(lat1)
      lat2 = radians(lat2)

      a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
      c = 2*asin(sqrt(a))

      return R * c

# Usage
lon1 = -103.548851
lat1 = 32.0004311
lon2 = -103.6041946
lat2 = 33.374939

print(haversine(lat1, lon1, lat2, lon2))

Solution 3

There is also a vectorized implementation, which allows to use 4 numpy arrays instead of scalar values for coordinates:

def distance(s_lat, s_lng, e_lat, e_lng):

   # approximate radius of earth in km
   R = 6373.0

   s_lat = s_lat*np.pi/180.0                      
   s_lng = np.deg2rad(s_lng)     
   e_lat = np.deg2rad(e_lat)                       
   e_lng = np.deg2rad(e_lng)  

   d = np.sin((e_lat - s_lat)/2)**2 + np.cos(s_lat)*np.cos(e_lat) * np.sin((e_lng - s_lng)/2)**2

   return 2 * R * np.arcsin(np.sqrt(d))

Solution 4

You can try the haversine package: https://pypi.org/project/haversine/

Example code:

from haversine import haversine
haversine((45.7597, 4.8422),(48.8567, 2.3508), unit='mi')
243.71209416020253

Solution 5

The bearing calculation is incorrect, you need to swap the inputs to atan2.

    bearing = atan2(sin(long2-long1)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(long2-long1))
    bearing = degrees(bearing)
    bearing = (bearing + 360) % 360

This will give you the correct bearing.

View more solutions

Comments

[eat]

Python haversine implementation can be found codecodex.com/wiki/…. However for short distance calculations very simple ways exists. Now, what is your maximum distance expected? Are you able to get your co-ordinates in some local cartesian co-ordinate system?

[Fábio Diniz]

Some implementations in python: - code.activestate.com/recipes/… - platoscave.net/blog/2009/oct/5/…

[eat]

@James Dyson: with distances like 15km, creat circle doesen't count anything. My suggestion: figure out first the solution with euclidean distances! That will give you a working solution and then later if your distances will be much much longer, then adjust your application. Thanks

[avitex]

Can I also find the Bearing from this equation?

[eat]

@James Dyson: If your above comment was aimed to me (and at to my earlier suggestion), the answer is surely (and quite 'trivially' as well). I may be able to give some example code, but it won't utilize trigonometry, rather geometry (so I'm unsure if it will help you at all. Are you familiar at all with the concept of vector? In your case positions and directions could be handled most straightforward manner with vectors).

[avitex]

@eat Sorry this message didn't get to you before...Yes it would be wonderful if you could give me some example code using the concept of vector. No I am not familiar with the concept vector yet.

[user102008]

atan2(sqrt(a), sqrt(1-a)) is the same as asin(sqrt(a))

[BasedRebel]

Could use math.radians() function instead of multiplying by pi/180 - same effect, but a bit more self-documenting.

[avitex]

Thanks soooo much, My distance can vary (horizontaly) between 15-1km

[avitex]

Correct me if I am wrong, but can't you just say: import math?

[Michael Dunn]

You can, but if you say import math then you have to specify math.pi, math.sin etc. With from math import * you get direct access to all the module contents. Check out "namespaces" in a python tutorial (such as docs.python.org/tutorial/modules.html)

[Eyal]

How come you use atan2(sqrt(a), sqrt(1-a)) instead of just asin(sqrt(a))? Is atan2 more accurate in this case?

[Michael Dunn]

No reason. I've tested this and there's no difference. I've changed the script to your simpler version.

[fangsterr]

Took me a while to figure out what the km = 6367 * c line meant. 6367 km is the radius of the Earth. If you want to do the same formula with miles, you should just multiply by 3956 instead.

[Dmitriy]

should be float division to cover really rare corner case of dlat|dlon being integers: a = sin(dlat/2.)**2 + cos(lat1) * cos(lat2) * sin(dlon/2.)**2

[Michael Dunn]

Good point, I hadn't thought of that. But in this case it's OK, since radians returns a float: radians(degrees(1)) gives 1.0.

[Gocht]

How can this be used in a Django's ORM query?

[Holger Bille]

I think it's just: Bearing = Bearing % 360

[ekhumoro]

If the mean Earth radius is defined as 6371 km, then that is equivalent to 3959 miles, not 3956 miles. See Global average radii for various ways to calculate these values.

[Alex van Es]

This one is much more accurate then the examples above!

[Chris Watts]

Altitude makes a difference to the radius by the way. If you're at an altitude of more than 300 metres I'd say it makes a significant difference, so factor it it.

[AesculusMaximus]

whats this returning? The bearing or the distance?

[Chris Owens]

@AesculusMaximus if using km for the radius of the earth it will return in kilometres distances between the two sets of coords. e.g. 0.062222488302154946 is 62.22 meters.

[ldmtwo]

This doesn't address the variation of R. 6356.752 km at the poles to 6378.137 km at the equator

[arilwan]

I am actually struggling to understand how these equations were derived as I am reading a paper. You have given me a pointer: haversine formula my first time to hear this, thank you.

[basilisk]

this method gives other results than all other methods above!

[Tejas Kale]

Does that error really matter to for your application? cs.nyu.edu/visual/home/proj/tiger/gisfaq.html

[Michael Dunn]

This would be a nice time to use augmented assignment too : Bearing %= 360

[Willem Hendriks]

great circle is usually implemented with haversine. So they could be the same with different radius.

[Anmol Deep]

Does haversine library have a function to calculate bearing?