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.
"""
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 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
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 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 getspeed's 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
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 kilometers per hour, which is passed to the decorator, which
converts the speed to the units specified.