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

Links | Site |
---|---|

GeographicLib | geographiclib.sourceforge.io |

Geodesics on an ellipsoid | en.wikipedia.org |

Haversine formula | en.wikipedia.org |

[How to co-locate a MODIS granule with CALIOP and CloudSat using python ? | |

](https://en.moonbooks.org/Jupyter/how_to_colocate_a_modis_granule_with_caliop_and_cloudsat/) | en.moonbooks.org |