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 "--------------------" 11 Answers
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 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)) 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* 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)) You can try the haversine package:
Example code:
from haversine import haversine haversine((45.7597, 4.8422),(48.8567, 2.3508), unit='mi') 243.71209416020253 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.
1Here'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 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.
2The 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 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 import statement to make the implementation.
Here are two functions to calculate distance and bearing, which are based on the code in previous messages and (added tuple type for geographical points in lat, lon format for both functions for clarity). I tested both functions and they seem to work right.
#coding: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) Considering that your goal is to measure the distance between two points (represented by geographic coordinates), will leave three options below:
Option 1
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 imports math.radians(x) and use them as follows
#Import radians from 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 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 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 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 in km, change `miles` to `km` [Out]: 536.997990696 
