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

Posted on

Solving problem is about exposing yourself to as many situations as possible like Haversine Formula in Python (Bearing and Distance between two GPS points) and practice these strategies over and over. With time, it becomes second nature and a natural way you approach any problems in general. Big or small, always start with a plan, use other strategies mentioned here till you are confident and ready to code the solution.
In this post, my aim is to share an overview the topic about Haversine Formula in Python (Bearing and Distance between two GPS points), which can be followed any time. Take easy to follow this discuss.

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

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 "--------------------"
Asked By: avitex

||

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 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
    return c * r
Answered By: Michael Dunn

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))
Answered By: Clay

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))
Answered By: Sergey Malyutin

Answer #4:

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.

Answered By: Jon Anderson

Answer #5:

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

Answer #6:

You can try the following:

from haversine import haversine
haversine((45.7597, 4.8422),(48.8567, 2.3508), unit='mi')
243.71209416020253
Answered By: Vamshi G

Answer #7:

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.

Answered By: OBu

Answer #8:

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)in(lat1)*cos(lat2)*cos(lon2-lon1))
bearing = degrees(bearing)
bearing = (bearing + 360) % 360
Answered By: gisdude

Leave a Reply

Your email address will not be published.