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 determines the great-circle
distance between two points on a sphere given their longitudes and
latitudes^{1}. The
formulaic representation of the Haversine formula is:

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})\): Geographic coordinate pair 1
- \((\phi_{2}, \lambda_{2})\): 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^{2}. 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)^{3}.
Due to this variation, the haversine distance calculation will always contain
some error, but for non-antipodal coordinate pairs provides a good
approximation. In the examples that follow, we specify a constant Earth radius
of **6367km**. Adjust as necessary for 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.

```
"""
Haversine formula implementation. For more information, refer to:
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
Radius of the Earth (est).
Returns
-------
float
Distance in kilometers between geoloc1 and geoloc2.
See Also
--------
.. [1] 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 Chicago to New York City, 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 resembles:

```
{
"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
Radius of the Earth (est).
Returns
-------
float
Distance in kilometers between geoloc1 and geoloc2.
See Also
--------
.. [1] 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 this 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 `getspeed`

‘s estimate, which isn’t bad considering our 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 specifies 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
declaration unmodified and wrapping the result returned by `legacyfunc`

with
logic to handle the change of units conversion would be the better choice.

We can implement 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
Radius of the Earth (est).
Returns
-------
float
Distance in kilometers between geoloc1 and geoloc2.
See Also
--------
.. [1] 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 miles per hour, which is then passed to the decorator, which
converts the speed to the units specified.

### Conclusion

The haversine formula can be applied in many cases, and as demonstrated in
this post, provides an acceptable method for computing distances between pairs
of geographic coordinates. As mentioned before, the constant Earth radius
assumption does not exactly model reality, but can still provide relatively
accurate distance measurements for locations not nearly antipodal.

We also introduced parameterized function decorators, explained their
motivation and demonstrated a real-world use case. Decorators have found wide
use in Python, and gaining a deeper understanding of their use cases will only
make you a better developer. Until next time, happy coding!