For many types of analysis, it is necessary to compute the distance between a pair of geographic coordinates, or the distance between multiple sets of coordinates and a central reference point. This can be readily accomplished by making use of the Haversine formula, which can be used to determine the great-circle distance between two points on a sphere given their longitudes and latitudes. The formulaic representation of the Haversine formula is

$$d = 2r\arcsin \left({\sqrt {\sin^{2}\left({\frac{\varphi_{2}- \varphi_{1}}{2}}\right)+\cos(\varphi_{1})\cos(\varphi_{2})\sin^{2} \left({\frac{\lambda_{2}-\lambda_{1}}{2}}\right)}}\right)$$

Where:

• $$\phi_{1}$$: Latitude of point 1
• $$\phi_{2}$$: Latitude of point 2
• $$\lambda_{1}$$: Longitude of point 1
• $$\lambda_{2}$$: Longitude of point 2
• $$r$$: The radius of the Earth

Using this notation:

• $$(\phi_{1}, \lambda_{1})$$ represents geographic coordinate pair 1
• $$(\phi_{2}, \lambda_{2})$$ represents geographic coordinate pair 2

Geographic coordinates can be specified in many different formats, but one common representation, and the one we’ll focus on for the remainder of the article is the decimal degrees format (e.g., {41.8781N, 87.6298W} for Chicago; {19.4326N, 99.1332W} for Mexico City, etc.).

Computer applications typically require the arguments of trigonometric functions to be expressed in radians. To convert decimal degrees to radians, multiply the number of degrees by $$\pi/180$$.

The Earth is not a perfect sphere, and as a result, the radius varies as a function of distance from the equator (the radius is 6357km at the poles and 6378km at the equator. Due to this variation, the haversine distance calculation will always contain some error, but for non-antipodal coordinate pairs still provides a good approximation. In the examples that follow, we specify a constant Earth radius of 6367km. Adjust as required to suit your purposes.

Typically, inverse trigonometric functions return results expressed in radians. This is the case for the Python Standard Library’s inverse trigonometric functions in the math library, which can be verified from the docstring for arcsin:

>>> import math
>>> help(math.asin)
Help on built-in function asin in module math:

asin(...)
asin(x)

Return the arc sine (measured in radians) of x.


To express the result in decimal degrees, multiply the number of radians returned by the inverse trigonometric function by $$180/\pi=57.295780$$ degrees/radian.

Next we create a function that computes the haversine formula and returns the distance between two geographic coordinate pairs specified in decimal degrees.

"""

http://www.movable-type.co.uk/scripts/gis-faq-5.1.html
"""
import math

def haverdist(geoloc1, geoloc2, R=6367):
"""
Compute distance between geographic coordinate pairs.

Parameters
----------
geoloc1: tuple or list
(lat1, lon1) of first geolocation.
geoloc2: tuple or list
(lat2, lon2) of second geolocation.
R: float

Returns
-------
float
Distance in kilometers between geoloc1 and geoloc2.

--------
..  http://www.movable-type.co.uk/scripts/gis-faq-5.1.html
"""
# Convert degress to radians then compute differences.
rlat1, rlon1 = [i * math.pi / 180 for i in geoloc1]
rlat2, rlon2 = [i * math.pi / 180 for i in geoloc2]
drlat, drlon = (rlat2 - rlat1), (rlon2 - rlon1)

init = (math.sin(drlat / 2.))**2 + (math.cos(rlat1)) * \
(math.cos(rlat2)) * (math.sin(drlon /2.))**2

# The intermediate result init is the great circle distance
# in radians. The return value will be in the same units as R.
# min protects against possible roundoff errors that could
# sabotage computation of the arcsine if the two points are
# very nearly antipodal, e.g., on opposide sides of the Earth.
return(2.0 * R * math.asin(min(1., math.sqrt(init))))


To test haverdist, we compute the distance between New York City and Chicago. For Chicago, use latitude=41.881832, and longitude=-87.623177 and New York City latitude=40.730610 and longitude=-73.935242. It’s now only a matter of passing the coordinates as tuples to the function and verifying the result:

>>> loc1 = (41.881832, -87.623177) # Chicago
>>> loc2 = (40.730610, -73.935242) # NYC
>>> haverdist(geoloc1=loc1, geoloc2=loc2, R=6367)
1148.514752233814


Thus, haverdist indicates that the distance between Chicago, Illinois and New York City is ~1149km.

If we compare this to Google Maps GeoDistance API using the same set of coordinate pairs, we get 1,296km. The discrepancy is due to the nature of the measurement: The distance calculated with haverdist computes the shortest distance between the two pairs of geographic coordinates, or the “air travel distance”. The Google Maps API calculates the driving distance, which is always greater than or equal to the air travel distance.

We can validate our measurement using distancefromto.net. Plugging in the distance from Chicago to New York, we see that the air travel distance is 1,149km (714 miles), which agrees with the distance computed using the haversine formula.

## Application

We can use the ability to calculate the distance between two or more geographic coordinate pairs for a variety of purposes. To demonstrate, we use the haverdist function to compute the velocity of the International Space Station with respect to the surface of the Earth, by taking two snapshots of the ISS’s Earth-relative position separated by a known time differential. The velocity is then the distance between the two sets of coordinates (from haverdist) divided by the time differential.

In order to obtain the position of the International Space Station relative to the Earth’s surface at any given moment, we leverage the Open Notify API, querying the following URL (which takes no inputs):

http://api.open-notify.org/iss-now.json


A successful response will look similar to the following:

{
"message": "success",
"timestamp": UNIX_TIME_STAMP,
"iss_position": {
"latitude": CURRENT_LATITUDE,
"longitude": CURRENT_LONGITUDE
}
}


According to the Open Notify website, polling should occur at 5 second intervals. We’ll use the same duration as the time differential component in the velocity calculation.

In what follows, we define two additional functions: getiss which returns the Open Notify json response with timestamped latitude and longitude of the ISS, and getspeed, which takes two coordinate pairs and the units in which the computed speed should be returned:

"""
Compute the speed of ISS relative to the surface of the Earth.
Distance calculation relies on the haversine function, which
estimates the radius of the Earth as a constant scalar.
"""
import datetime
import math
import os
import sys
import time
import requests

def getiss(proxy=None):
"""
Get timestamped geo-coordinates of International Space Station.

Parameters
----------
proxy: dict, optional
Proxy details if required to make web requests.
Defaults to None.

Returns
-------
dict
Dictionary with keys "latitude", "longitude" and
"timestamp" indicating time and position of ISS.
"""
dpos = dict()
url  = "http://api.open-notify.org/iss-now.json"
resp = requests.get(url, proxies=proxy).json()
if resp["message"]=="success":
dpos["timestamp"] = resp["timestamp"]
dpos["latitude"]  = float(resp["iss_position"]["latitude"])
dpos["longitude"] = float(resp["iss_position"]["longitude"])
else:
raise RuntimeError("Unable to access Open Notify API")
return(dpos)

def haverdist(geoloc1, geoloc2, R=6367):
"""
Compute distance between geographic coordinate pairs in kilometers.

Parameters
----------
geoloc1: tuple or list
(lat1, lon1) of first geolocation.
geoloc2: tuple or list
(lat2, lon2) of second geolocation.
R: float

Returns
-------
float
Distance in kilometers between geoloc1 and geoloc2.

--------
..  http://www.movable-type.co.uk/scripts/gis-faq-5.1.html
"""
# Convert degress to radians then compute differences.
rlat1, rlon1 = [i * math.pi / 180 for i in geoloc1]
rlat2, rlon2 = [i * math.pi / 180 for i in geoloc2]
drlat, drlon = (rlat2 - rlat1), (rlon2 - rlon1)

init = (math.sin(drlat / 2.))**2 + (math.cos(rlat1)) * \
(math.cos(rlat2)) * (math.sin(drlon /2.))**2

# The intermediate result init is the great circle distance
# in radians. The return value will be in the same units as R.
# min protects against possible roundoff errors that could
# sabotage computation of the arcsine if the two points are
# very nearly antipodal, e.g., on opposide sides of the Earth
return(2.0 * R * math.asin(min(1., math.sqrt(init))))

def getspeed(dloc1, dloc2):
"""
Compute speed of ISS relative to Earth's surface using
a pair of coordinates retrieved from getiss.

Parameters
----------
dloc1: dict
Dictionary with keys "latitude", "longitude" "timestamp"
associated with the first positional snapshot.
dloc2: dict
Dictionary with keys "latitude", "longitude" "timestamp"
associated with the second positional snapshot.

Returns
-------
float
Scalar value representing the average speed of the International
Space Station relative to the Earth going from dloc1 to
dloc2.
"""
# Convert unix epochs to timestamp datetime objects.
ts1   = datetime.datetime.fromtimestamp(dloc1['timestamp'])
ts2   = datetime.datetime.fromtimestamp(dloc2['timestamp'])
secs  = abs((ts2-ts1).total_seconds())
loc1  = (dloc1["latitude"], dloc1["longitude"])
loc2  = (dloc2["latitude"], dloc2["longitude"])
dist  = haverdist(geoloc1=loc1, geoloc2=loc2)
vinit = (dist / secs) # kilometers per second
return(vinit * 3600)

# Typically, the following would be collected under a main
# function and called from the command line, but for ease of
# demonstration, we present the logic immediately following
# the function declarations.
ts_init = datetime.datetime.now().strftime("%c")
print(f"[{ts_init}] Retrieving ISS geographic coordinate snapshots...")
dpos1 = getiss()
time.sleep(5)
dpos2 = getiss()
speed = getspeed(dloc1=dpos1, dloc2=dpos2)
ts_final = datetime.datetime.now().strftime("%c")
print(f"[{ts_final}] ISS speed relative to Earth's surface: {speed:.2f}km/h")


Running the above produces the following:

[Sat Dec 29 12:40:30 2018] Retrieving ISS geographic coordinate snapshots...
[Sat Dec 29 12:40:35 2018] ISS speed relative to Earth's surface: 26253.60km/h


According to our calculations, the average orbital speed of the ISS relative to the Earth’s surface for the targeted time interval exceeded 26,000km/h (~16,313.23mph). Wikipedia’s ISS page puts the average orbital speed at 27,600km/h (17,100 mph), a difference of only 5% when compared with getspeeds estimate, which isn’t bad considering the constant Earth radius assumption.

## A Note on Parameterized Decorators

The official definition of a decorator is a function that takes another function and extends the behavior of the latter function without explicitly modifying it. I never really understood the purpose of decorators until coming to the realization that they can be used as a way to update the units of a value returned by a function without needing to alter the nature of the function of interest in anyway.

To illustrate, we refer back to getspeed declared in the previous code block. This function returns a scalar representing the average speed of the ISS over a given time differential. The result is returned in kilometers per hour, with no parameter available to apply a change of units. To incorporate this functionality, we can either (1) add a units parameter to change the units of the calculated speed prior to returning to caller, or (2) implement a decorator which converts the units of the calculated value independent of getspeed.

In this case, the first option seems like an obvious choice, but instead of getspeed imagine a different function (we’ll call it legacyfunc) that we didn’t author, which has been in production for a very long time, which is very difficult step through due to overly-complicated control structures, has lots of strange optional parameters and many more lines of code than getspeed, was responsible for returning the value of interest, and it is this value that requires a change of units. In this case, leaving the original callable unmodified and wrapping the result returned by legacyfunc with logic to handle the units conversion would be the better choice.

We can create a function to handle the change of units conversion. This will result in a parameterized decorator, the parameter indicating which units to return the orbital speed. For this implementation, the only options will be kilometers per hour or miles per hour, but the decorator can be extended to facilitate any number of additional distance or time conversions.

"""
Demonstration of parameterized decorator to facilitate change of
units conversion.
"""
import functools

def units(spec):
"""
Specify the units to represent orbital speed. spec can be either
"kmph" (kilometer per hour) or "mph" (miles per hour). Defaults
to "kmph".
"""
def decorator(function):
@functools.wraps(function)
def wrapper(*args, **kwargs):
resinit = function(*args, **kwargs)
# resinit is speed represented in kilometers per hour.
if spec=="mph":
result = resinit * 0.62137119
else:
result = res_init
return(result)
return(wrapper)
return(decorator)


Once units has been declared, we modify the declaration of getspeed by referencing the units decorator as follows:

import functools

@units("mph")
def getspeed(dloc1, dloc2):
"""
Compute speed of ISS relative to Earth's surface using
a pair of coordinates retrieved from getiss.

Parameters
----------
dloc1: dict
Dictionary with keys "latitude", "longitude" "timestamp"
associated with the first positional snapshot.
dloc2: dict
Dictionary with keys "latitude", "longitude" "timestamp"
associated with the second positional snapshot.

Returns
-------
float
Scalar value representing the average speed of the International
Space Station relative to the Earth going from dloc1 to
dloc2.
"""
# Convert unix epochs to timestamp datetime objects.
ts1   = datetime.datetime.fromtimestamp(dloc1['timestamp'])
ts2   = datetime.datetime.fromtimestamp(dloc2['timestamp'])
secs  = abs((ts2-ts1).total_seconds())
loc1  = (dloc1["latitude"], dloc1["longitude"])
loc2  = (dloc2["latitude"], dloc2["longitude"])
dist  = haverdist(geoloc1=loc1, geoloc2=loc2)
vinit = (dist / secs) # kilometers per second
return(vinit * 3600)


With this change, the scalar representing kilometers per hour returned by getspeed will be converted to miles per hour. Here’s the full implementation repeated, the differences being the units decoarator and decorated getspeed function:

import datetime
import functools
import math
import os
import sys
import time
import requests

def units(spec):
"""
Specify the units to represent orbital speed. spec can be either
"kmph" (kilometer per hour) or "mph" (miles per hour). Defaults
to "kmph".
"""
def decorator(function):
@functools.wraps(function)
def wrapper(*args, **kwargs):
resinit = function(*args, **kwargs)
# resinit is speed represented in kilometers per hour
if spec=="mph":
result = resinit * 0.62137119
else: result = res_init
return(result)
return(wrapper)
return(decorator)

def getiss(proxy=None):
"""
Get timestamped geo-coordinates of International Space Station.

Parameters
----------
proxy: dict, optional
Proxy details if required to make web requests.
Defaults to None.

Returns
-------
dict
Dictionary with keys "latitude", "longitude" and
"timestamp" indicating time and position of ISS.
"""
dpos = dict()
url  = "http://api.open-notify.org/iss-now.json"
resp = requests.get(url, proxies=proxy).json()
if resp["message"]=="success":
dpos["timestamp"] = resp["timestamp"]
dpos["latitude"]  = float(resp["iss_position"]["latitude"])
dpos["longitude"] = float(resp["iss_position"]["longitude"])
else:
raise RuntimeError("Unable to access Open Notify API")
return(dpos)

def haverdist(geoloc1, geoloc2, R=6367):
"""
Compute distance between geographic coordinate pairs in kilometers.

Parameters
----------
geoloc1: tuple or list
(lat1, lon1) of first geolocation.
geoloc2: tuple or list
(lat2, lon2) of second geolocation.
R: float

Returns
-------
float
Distance in kilometers between geoloc1 and geoloc2.

--------
..  http://www.movable-type.co.uk/scripts/gis-faq-5.1.html
"""
# Convert degress to radians then compute differences.
rlat1, rlon1 = [i * math.pi / 180 for i in geoloc1]
rlat2, rlon2 = [i * math.pi / 180 for i in geoloc2]
drlat, drlon = (rlat2 - rlat1), (rlon2 - rlon1)

init = (math.sin(drlat / 2.))**2 + (math.cos(rlat1)) * \
(math.cos(rlat2)) * (math.sin(drlon /2.))**2

# The intermediate result init is the great circle distance
# in radians. The return value will be in the same units as R.
# min protects against possible roundoff errors that could
# sabotage computation of the arcsine if the two points are
# very nearly antipodal, e.g., on opposide sides of the Earth
return(2.0 * R * math.asin(min(1., math.sqrt(init))))

@units("mph")
def getspeed(dloc1, dloc2):
"""
Compute speed of ISS relative to Earth's surface using
a pair of coordinates retrieved from getiss.

Parameters
----------
dloc1: dict
Dictionary with keys "latitude", "longitude" "timestamp"
associated with the first positional snapshot.
dloc2: dict
Dictionary with keys "latitude", "longitude" "timestamp"
associated with the second positional snapshot.

Returns
-------
float
Scalar value representing the average speed of the International
Space Station relative to the Earth going from dloc1 to
dloc2.
"""
# Convert unix epochs to timestamp datetime objects.
ts1   = datetime.datetime.fromtimestamp(dloc1['timestamp'])
ts2   = datetime.datetime.fromtimestamp(dloc2['timestamp'])
secs  = abs((ts2-ts1).total_seconds())
loc1  = (dloc1["latitude"], dloc1["longitude"])
loc2  = (dloc2["latitude"], dloc2["longitude"])
dist  = haverdist(geoloc1=loc1, geoloc2=loc2)
vinit = (dist / secs) # kilometers per second
return(vinit * 3600)

# Note that as a result of decorating getspeed with the units,
# the ISS speed estimate will be in miles per hour instead of
# kilometers per hour.
ts_init = datetime.datetime.now().strftime("%c")
print(f"[{ts_init}] Retrieving ISS geographic coordinate snapshots...")
dpos1 = getiss()
time.sleep(5)
dpos2 = getiss()
speed = getspeed(dloc1=dpos1, dloc2=dpos2)
ts_final = datetime.datetime.now().strftime("%c")
print(f"[{ts_final}] ISS speed relative to Earth's surface: {speed:.2f}mph")


Running results in:

[Sat Dec 29 13:44:17 2018] Retrieving ISS geographic coordinate snapshots...
[Sat Dec 29 13:44:22 2018] ISS speed relative to Earth's surface: 16126.21mph


Once decorated, the call to getspeed is no different than before. The speed is returned in kilometers per hour, which is passed to the decorator, which converts the speed to the units specified.