### Parametric Thoughts Image from New Old Stock

# Calculate Distance Between GPS Points in Python

When working with GPS, it is sometimes helpful to calculate distances between points. But simple Euclidean distance doesn’t cut it since we have to deal with a sphere, or an oblate spheroid to be exact. So we have to take a look at geodesic distances.

# Haversine Distance Formula

There are various ways to handle this calculation problem. For example there is the Great-circle distance, which is the shortest distance between two points on the surface of a sphere. Another similar way to measure distances is by using the Haversine formula, which takes the equation

\begin{equation} a = hav(\Delta\varphi) + cos(\varphi_1) \cdot cos(\varphi_2) \cdot hav(\Delta\lambda) \end{equation}

with haversine function

\begin{equation} hav(\theta) = sin^{2}(\frac{\theta}{2}) \end{equation}

and takes this to calculate the geodesic distance

\begin{equation} \text{distance} = 2 \cdot R \cdot arctan(\sqrt{a}, \sqrt{1-a}) \end{equation}

where the latitude is $\varphi$, the longitude is denoted as $\lambda$ and $R$ corresponds to Earths mean radius in kilometers which is approximately 6371 km. We can take this formula now and translate it into Python:

import math

def haversine(coord1, coord2):
R = 6372800  # Earth radius in meters
lat1, lon1 = coord1
lat2, lon2 = coord2

a = math.sin(dphi/2)**2 + \
math.cos(phi1)*math.cos(phi2)*math.sin(dlambda/2)**2

return 2*R*math.atan2(math.sqrt(a), math.sqrt(1 - a))


Important to note is that we have to take the radians of the longitude and latitude values. We can take this function now and apply distances to different cities. Lets say we want to calculate the distances from London to some other cities.

london_coord = 51.5073219,  -0.1276474
cities = {
'berlin': (52.5170365,  13.3888599),
'vienna': (48.2083537,  16.3725042),
'sydney': (-33.8548157, 151.2164539),
}

for city, coord in cities.items():
distance = haversine(london_coord, coord)
print(city, distance)

berlin 930723.2019867424
vienna 1235650.1412429416
sydney 16997984.55171465


This already gives us some seemingly accurate result, but let’s compare it to another method.

You can also use geopy to measure distances. This package has many different methods for calculating distances, but it uses the Vincenty’s formulae as default, which is a more exact way to calculate distances on earth since it takes into account that the earth is, as previously mentioned, an oblate spheroid. The Vincenty’s formulae is well described in this article.

from geopy.distance import distance

for city, coord in cities.items():
d = distance(london_coord, coord).m
print(city, d)

berlin 933410.7641236288
vienna 1238804.7757673298
sydney 16988546.466908157


As you can see, there is a difference between the values, especially since we work with very large distances, which enhances the distortion of our spheroid-shaped Earth.

# Projections Using pyproj

There is also the pyproj Python package, which offers Python interfaces to PROJ.4. It is a great package to work with map projections, but in there you have also the pyproj.Geod class which offers various geodesic computations. To calculate the distance between two points we use the pyproj.Geod.inv function, which calculates an inverse transformation and returns forward and back azimuths and distance.

import pyproj

geod = pyproj.Geod(ellps='WGS84')

for city, coord in cities.items():
lat0, lon0 = london_coord
lat1, lon1 = coord

azimuth1, azimuth2, distance = geod.inv(lon0, lat0, lon1, lat1)
print(city, distance)
print('    azimuth', azimuth1, azimuth2)

berlin 933410.7641236288
azimuth 77.79312482066598 -91.53477000281634
vienna 1238804.7757673291
azimuth 100.7430617124229 -66.6018639905512
sydney 16988546.466908157
azimuth 60.33221400668487 -40.68498881273351