Calculate Distance Between GPS Points in Python

Table of Contents

• Haversine Distance Formula
• Projections Using pyproj

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
    
    phi1, phi2 = math.radians(lat1), math.radians(lat2) 
    dphi       = math.radians(lat2 - lat1)
    dlambda    = math.radians(lon2 - lon1)
    
    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),
    'madrid': (40.4167047,  -3.7035825) 
}

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

berlin 930723.2019867424
vienna 1235650.1412429416
sydney 16997984.55171465
madrid 1263769.8859593808

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

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
madrid 1263101.92391795
    azimuth -166.0130675331932 11.403752336198153

On a geographic sidenote, the forward azimuth is the direction which is defined as a horizontal angle measured clockwise from a north base line and a back azimuth is the opposite direction of the forward azimuth. You could use this information for example to sail the ocean if this is what you intend.