How to calculate the distance between two geographical points defined by latitudes and longitudes using Python?

Introduction

In this tutorial, we will be learning how to calculate the distance between two geographical points using python. This can be useful in many scenarios, such as finding the shortest route for a road trip or determining the closest city from your current location.

To begin with, it is important to understand that geographical coordinates are measured in latitude and longitude. Latitude represents the distance north or south of the equator, while longitude represents the distance east or west of the prime meridian.

Now let's dive into the steps for calculating distance between two geographical points using python:

Calculating Distance between Geographical Points using Python

Using geographiclib

To calculate the distance between two geographical points defined by latitude and longitude using Python, a first solution is to use the geographiclib library which provides a simple interface for performing geometric calculations on an ellipsoidal model of the earth.

To define the two geographical points for which we want to calculate the distance between. We will be using latitude and longitude coordinates for this task. Let's say we want to caclulate the distance between New-York city and Paris:

ny_lat = 40.7128
ny_lon = - 74.0060

paris_lat = 48.8566
paris_lon = 2.3522

Now that we have our points defined, we can use the Geodesic class from the geographiclib library to calculate the distance between them:

from geographiclib.geodesic import Geodesic

g = Geodesic.WGS84.Inverse(ny_lat, ny_lon, paris_lat, paris_lon)

print(g)

The output of this calculation will be a dictionary containing various values such as the distance in meters, kilometers, and miles:

{'lat1': 40.7128, 'lon1': -74.006, 'lat2': 48.8566, 'lon2': 2.3522, 'a12': 52.653277894384246, 's12': 5852935.2917667655, 'azi1': 53.72404460168419, 'azi2': 111.82608319619803}

To get the distance in meters, we can access the s12 key from the dictionary:

d = g['s12']

Then

print("Distance is {:.2f}m".format(d))

print("Distance is {:.2f}km".format(d/1000.0))

gives

Distance is 5852935.29m
Distance is 5852.94km

respectively.

Implementing the Haversine Formula in Python

An alternative approach is to personally implement the Haversine Formula in Python, which can be sourced from platforms like Wikipedia. This method allows for a more hands-on experience and a deeper understanding of the formula itself.

The Haversine Formula is a mathematical equation used to calculate the distance between two points on a sphere. It takes into account the curvature of Earth and provides accurate results for long distances. The formula uses the latitude and longitude coordinates of two points to calculate the distance between them:

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

ny_coord = (ny_lat, ny_lon)
paris_coord = (paris_lat, paris_lon)

d = haversine(ny_coord, paris_coord)

print("Distance is {:.2f}m".format(d))

print("Distance is {:.2f}km".format(d/1000.0))

gives

Distance is 5838890.10m
Distance is 5838.89km

respectively.

Distance between an array of coordinates

Sometimes, calculating the distance between two points is not enough. We often need to calculate the distance between a collection of coordinates stored in a numpy array. For such cases, we can leverage the power of Python's broadcasting approach. This allows us to efficiently perform calculations on arrays, enabling us to handle complex data with ease.

Example of application in remote sensing: How to co-locate a MODIS granule with CALIOP and CloudSat using python ?

import numpy as np

def broadcasting_based_distance(data1, data2):

    # The arrays data1 and data2 contain latitude and longitude values respectively, 
    # represented in two columns

    R = 6372.8  # Earth radius in kilometers

    data1 = np.deg2rad(data1)                     
    data2 = np.deg2rad(data2)

    lat1 = data1[:,0]                     
    lng1 = data1[:,1]

    lat2 = data2[:,0]                     
    lng2 = data2[:,1]

    diff_lat = lat1[:,None] - lat2
    diff_lng = lng1[:,None] - lng2
    d = np.sin(diff_lat/2)**2 + np.cos(lat1[:,None])*np.cos(lat2) * np.sin(diff_lng/2)**2

    return 2 * 6373 * np.arcsin(np.sqrt(d))


london_coord = np.array( [[ny_lat, ny_lon]] )
madrid_coord = np.array( [[paris_lat, paris_lon]] )

print(london_coord.shape)

D = broadcasting_based_distance(london_coord,madrid_coord)

d = D[0][0]

print(d)

print("Distance is {:.2f}km".format(d))

returns

Distance is 5839.07km

Plotting the distance between two geographical points

from matplotlib.pyplot import figure

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import matplotlib.patches as mpatches

fig = figure(num=None, figsize=(12, 10), dpi=100, edgecolor='k')

map_proj = ccrs.PlateCarree()

ny_lat = 40.7128
ny_lon = - 74.0060

paris_lat = 48.8566
paris_lon = 2.3522

ax = plt.axes(projection=map_proj)

ax.set_global()
ax.gridlines()

ax.coastlines(linewidth=0.5, color='k', resolution='50m')

text_lon = ny_lon
text_lat = ny_lat

text = "New-York Paris \n Distance \n {:.2f}km".format(d)

ax.text((text_lon+paris_lon)/2.0 - 30.0,
        text_lat - 20.0,
        text,
        fontsize=12,
        color='coral',
        weight='bold')

plt.plot([ny_lon, paris_lon], [ny_lat, paris_lat],
         color='coral', linewidth=2, marker='o',transform=ccrs.Geodetic())

plt.title('How to calculate the distance between two geographical points \n defined by latitudes and longitudes using Python?')

plt.savefig("distance_01.png", bbox_inches='tight', facecolor='white')

plt.show()

References