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 two GPS points.

I have researched on the haversine distance. Someone told me that I could also find the bearing using the same data.

---

Everything is working fine, but the bearing doesn't quite work right yet. The bearing outputs negative, but it 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 "--------------------"

Revelry

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? – Drogheda

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

@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 – Drogheda

Can I also find the Bearing from this equation? – Revelry

@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). – Drogheda

@Drogheda 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. – Revelry

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

Answer 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

Manofwar

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

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

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) – Manofwar

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

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

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. – Julide

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 – Useful

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. – Manofwar

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. – Salep

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. – Feminize

whats this returning? The bearing or the distance? – Patino

@Patino 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. – Thaumaturge

Answer 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))

Dominant

This one is much more accurate then the examples above! – Annotate

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

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

@AlexvanEs Adding precision with bad accuracy has no added benefit. In fact it's the opposite since it gives a false sense of good accuracy https://mcmap.net/q/111695/-what-is-the-difference-between-39-precision-39-and-39-accuracy-39-closed – Tranquil

Re "these functions will be off": Can you quantify that? 10%? 1%? 0.000001%? – Bulletproof

Answer 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))

Dropsy

Answer 4

You can try the haversine package:

Example code:

from haversine import haversine

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

Output:

243.71209416020253

Delcine

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

Does haversine library have a function to calculate bearing? – Yuille

Answer 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.

Ruben

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. – Lockup

Answer 6

Here's a NumPy vectorized implementation of the Haversine Formula given by @Michael Dunn, gives a 10-50 times improvement over large vectors.

from numpy import radians, cos, sin, arcsin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    """

    # Convert decimal degrees to radians:
    lon1 = np.radians(lon1.values)
    lat1 = np.radians(lat1.values)
    lon2 = np.radians(lon2.values)
    lat2 = np.radians(lat2.values)

    # Implementing the haversine formula:
    dlon = np.subtract(lon2, lon1)
    dlat = np.subtract(lat2, lat1)

    a = np.add(np.power(np.sin(np.divide(dlat, 2)), 2),
               np.multiply(np.cos(lat1),
                           np.multiply(np.cos(lat2),
                               np.power(np.sin(np.divide(dlon, 2)), 2))))

    c = np.multiply(2, np.arcsin(np.sqrt(a)))
    r = 6371

    return c*r

Simonesimoneau

This code doesn't work; you import all the functions from NumPy and then you use np all over your code. – Boyse

Answer 7

Considering that your goal is to measure the distance between two points (represented by geographic coordinates), will leave three options below:

1. Haversine formula

2. Using GeoPy geodesic distance

3. Using GeoPy great-circle distance

---

Option 1

The haversine formula will do the work. However, it is important to note that by doing that one is approximating the Earth as a sphere, and that has an error (see this answer) - as Earth is not a sphere.

In order to use the haversine formula, first of all, one needs to define the radius of the Earth. This, in itself, may lead to some controversy. Considering the following three sources

• NASA's Goddard Space Flight Center: 6371 km

• Wikipedia: 6371 km (3958.8 mi)

• Google - 6371 km

I'll be using the value 6371 km as a reference to the radius of the Earth.

# Radius of the Earth
r = 6371.0

We will be leveraging math module.

After the radius, one moves to the coordinates, and one starts by converting the coordinated into radians, in order to use math's trigonometric functions. For that one, it imports math.radians(x) and use them as follows:

# Import radians from the 'math' module
from math import radians

# Latitude and longitude for the first point (let's consider 40.000º and 21.000º)
lat1 = radians(40.000)
lon1 = radians(21.000)

# Latitude and longitude for the second point (let's consider 30.000º and 25.000º)
lat2 = radians(30.000)
lon2 = radians(25.000)

Now one is ready to apply the haversine formula. First, one subtracts the longitude of point 1 to the longitude of point 2

dlon = lon2 - lon1
dlat = lat2 - lat1

Then, and for here there are a couple of trigonometric functions that one is going to use, more specifically, math.sin(), math.cos(), and math.atan2(). We will also be using math.sqrt()

# Import sin, cos, atan2, and sqrt from the 'math' module
from math import sin, cos, atan2, sqrt

a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
d = r * c

Then one gets the distance by printing d.

As it may help, let's gather everything in a function (inspired by Michael Dunn's answer)

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

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great-circle distance (in km) between two points
    using their longitude and latitude (in degrees).
    """
    # Radius of the Earth
    r = 6371.0

    # Convert degrees to radians
    # First point
    lat1 = radians(lat1)
    lon1 = radians(lon1)

    # Second point
    lat2 = radians(lat2)
    lon2 = radians(lon2)

    # Haversine formula
    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))
    return r * c

---

Option 2

One is going to use GeoPy's distance, more specifically, the geodesic.

We can obtain the results both on km, or miles (Source)

# Import Geopy's distance
from geopy import distance

wellington = (-41.32, 174.81)
salamanca = (40.96, -5.50)
print(distance.distance(wellington, salamanca).km) # If one wants it in miles, change `km` to `miles`

[Out]: 19959.6792674

---

Option 3

One is going to use GeoPy's distance, more specifically, the great-circle.

We can obtain the results both on km, or miles (Source)

# Import Geopy's distance
from geopy import distance

newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)

print(distance.great_circle(newport_ri, cleveland_oh).miles) # If one wants it in km, change `miles` to `km`

[Out]: 536.997990696

---

Notes:

• As the great-circle distance is often calculated using the Haversine formula (as Willem Hendriks noted), Option 1 and 3 are similar, but use a different radius.

  • GeoPy's Great-circle distance uses a spherical model of the earth, using the mean earth radius as defined by the International Union of Geodesy and Geophysics, 6371.0087714150598 kilometers approx. 6371.009 km (for WGS-84), resulting in an error of up to about 0.5% [Source].

Gerdagerdeen

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

Answer 8

You can solve the negative bearing problem by adding 360°. Unfortunately, this might result in bearings larger than 360° for positive bearings. This is a good candidate for the modulo operator, so all in all you should add the line

Bearing = (Bearing + 360) % 360

at the end of your method.

Illume

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

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

Answer 9

The Y in atan2 is, by default, the first parameter. Here is the documentation. You will need to switch your inputs to get the correct bearing angle.

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

Aberdare

This looks like a bogus answer. Is it? What is "in" supposed to do? – Bulletproof

@PeterMortensen Yes it is/was... Looks like there has been a typo in there for the past 7 years... I don't remember it being wrong last time I looked though... ? Thank you for pointing it out. I edited it. – Aberdare

Answer 10

Refer to Difference between Vincenty and great-circle distance calculations.

This actually gives two ways of getting distance. They are haversine and Vincentys. From my research, I came to know that Vincentys is relatively accurate. Also use an import statement to make the implementation.

Saboteur

Answer 11

Here are two functions to calculate distance and bearing, which are based on the code in previous messages and Compass bearing between two points in Python (I added a tuple type for geographical points in latitude, longitude format for both functions for clarity). I tested both functions, and they seemed to work right.

# Encoding: UTF-8
from math import radians, cos, sin, asin, sqrt, atan2, degrees

def haversine(pointA, pointB):

    if (type(pointA) != tuple) or (type(pointB) != tuple):
        raise TypeError("Only tuples are supported as arguments")

    lat1 = pointA[0]
    lon1 = pointA[1]

    lat2 = pointB[0]
    lon2 = pointB[1]

    # Convert decimal degrees to radians
    lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2])

    # 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
    return c * r

def initial_bearing(pointA, pointB):

    if (type(pointA) != tuple) or (type(pointB) != tuple):
        raise TypeError("Only tuples are supported as arguments")

    lat1 = radians(pointA[0])
    lat2 = radians(pointB[0])

    diffLong = radians(pointB[1] - pointA[1])

    x = sin(diffLong) * cos(lat2)
    y = cos(lat1) * sin(lat2) -
        (sin(lat1) * cos(lat2) * cos(diffLong))

    initial_bearing = atan2(x, y)

    # Now we have the initial bearing but math.atan2 return values
    # from -180° to + 180° which is not what we want for a compass bearing
    # The solution is to normalize the initial bearing as shown below
    initial_bearing = degrees(initial_bearing)
    compass_bearing = (initial_bearing + 360) % 360

    return compass_bearing

pA = (46.2038, 6.1530)
pB = (46.449, 30.690)

print haversine(pA, pB)

print initial_bearing(pA, pB)

Tashia

this method gives other results than all other methods above! – Kilo

Answer 12

You can use the below implementation in Python

import math

def haversine_distance(lat1, lon1, lat2, lon2, unit='K'):
    r = 6371  # radius of the earth in kilometers
    if unit == 'M':
        r = 3960  # radius of the earth in miles
    dLat = math.radians(lat2 - lat1)
    dLon = math.radians(lon2 - lon1)
    a = math.sin(dLat / 2) * math.sin(dLat / 2) + \
        math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * \
    math.sin(dLon / 2) * math.sin(dLon / 2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    distance = r * c
    return distance

You can read more about it at Haversine Formula

Urbain