Calculating Great Circle Distances in Python
The shortest distance between two locations on the surface of Earth (or any planet) is known as the Great Circle Distance. Except for relatively short distances these cannot be measured on a map due to the distortion and flattening necessary to represent a sphere on a flat plane. However, the calculation of the Great Circle Distance from two coordinates is simple, although I suspect generations of midshipmen might not have agreed. In this post I will write a short program in Python to calculate the distance from London to various other cities round the world.
As always, Wikipedia provides a comprehensive reference on the subject of calculating Great Circle distances. The principle is very simple, if you understand radians.
Everyone is familiar with using degrees to measure angles - we all gain an understanding of them when we are given our first school protractor. Slightly less well understood is the radian. (If you already know all about radians skip to the next paragraph!) Imagine a circle (for the purposes of this project, imagine the circle is the flat part of a hemisphere, a sphere cut in half) and now imagine drawing a line from the centre to the edge. The length of that line is of course the radius. Now draw round the circumference of the circle by the same distance as the radius, finally drawing another line back to the centre of the circle. The angle between the two lines meeting at the centre is 1 radian. The size of 1 radian in degrees is 180/π which is an irrational number but 57.29577951° to 8dp.
If you skipped my scholarly discourse on radians, welcome back. If you think about it for a moment, it should become clear that if the angle between two points on the surface of our little blue planet is 1 radian, then the great circle distance between them is the same as the radius of the planet. If the angle is 2 radians, the distance is twice the radius; if the angle is 0.5 radians the distance is half the radius and so in. In short, the Great Circle Distance between two points in the angle between them multiplied by the radius. The units of distance are irrelevant as long as they are the same (don't mix miles and kilometres for example) but we must use radians rather than degrees. This is the reason we are using radians - they reduce the problem to one simple multiplication.
We're getting ahead of ourselves a bit here, as we first need to calculate the angle between two points from their latitudes and longitudes. I'll "borrow" the formula for this from Wikipedia, but first let's define some terms.
Φ1 and λ1 - latitude and longitude in radians of first location
Φ2 and λ2 - latitude and longitude in radians of second location
ΔΦ and Δλ - absolute differences between latitudes and longitudes
Δσ - angle between locations
Don't worry, I won't be using Greek in the code. Now for the formula - note use of dot notation for multiplication.
So basically we need to apply this formula to a pair of coordinates, then multiply the result by the radius of the Earth.
The Code
The code consists of the following two files which you can grab from the Github repository.
greatcircle.py
main.py
This is greatcircle.py.
import math
DEGREES_IN_RADIAN = 57.29577951
MEAN_EARTH_RADIUS_KM = 6371
KILOMETRES_IN_MILE = 1.60934
class GreatCircle(object):
"""
This class has attributes for the names and locations of a pair of cities,
and for their distances.
After the city attributes are set, call calculate to set the distance attributes.
Validation is carried out and the valid attributes set. This should be checked before using
calculated attributes.
"""
def __init__(self):
"""
Create a set of attributes with default values.
"""
self.name1 = None
self.latitude1_degrees = 0
self.longitude1_degrees = 0
self.latitude1_radians = 0
self.longitude1_radians = 0
self.name2 = None
self.latitude2_degrees = 0
self.longitude2_degrees = 0
self.latitude2_radians = 0
self.longitude2_radians = 0
self.central_angle_radians = 0
self.central_angle_degrees = 0
self.distance_kilometres = 0
self.distance_miles = 0
self.valid = False
def calculate(self):
"""
Central method to set calculated attributes, which it
does by calling other private functions.
"""
self.__validate_degrees()
if self.valid:
self.__calculate_radians()
self.__calculate_central_angle()
self.__calculate_distance()
def __validate_degrees(self):
"""
Check latitudes and longitudes are within valid ranges,
setting the valid attribute accordingly.
"""
self.valid = True
if self.latitude1_degrees < -90.0 or self.latitude1_degrees > 90.0:
self.valid = False
if self.longitude1_degrees < -180.0 or self.longitude1_degrees > 180.0:
self.valid = False
if self.latitude2_degrees < -90.0 or self.latitude2_degrees > 90.0:
self.valid = False
if self.longitude2_degrees < -180.0 or self.longitude2_degrees > 180.0:
self.valid = False
def __calculate_radians(self):
"""
Calculate radians from degrees by dividing by constant.
"""
self.latitude1_radians = self.latitude1_degrees / DEGREES_IN_RADIAN
self.longitude1_radians = self.longitude1_degrees / DEGREES_IN_RADIAN
self.latitude2_radians = self.latitude2_degrees / DEGREES_IN_RADIAN
self.longitude2_radians = self.longitude2_degrees / DEGREES_IN_RADIAN
def __calculate_central_angle(self):
"""
Slightly complex formula for calculating the central angle
between two points on the surface of a sphere.
"""
if self.longitude1_radians > self.longitude2_radians:
longitudes_abs_diff = self.longitude1_radians - self.longitude2_radians
else:
longitudes_abs_diff = self.longitude2_radians - self.longitude1_radians
self.central_angle_radians = math.acos( math.sin(self.latitude1_radians)
* math.sin(self.latitude2_radians)
+ math.cos(self.latitude1_radians)
* math.cos(self.latitude2_radians)
* math.cos(longitudes_abs_diff))
self.central_angle_degrees = self.central_angle_radians * DEGREES_IN_RADIAN
def __calculate_distance(self):
"""
Because we are using radians, this is a simple formula multiplying the radius
by the angle, the actual units used being irrelevant.
Also the distance in miles is calculated from kilometres.
"""
self.distance_kilometres = MEAN_EARTH_RADIUS_KM * self.central_angle_radians
self.distance_miles = self.distance_kilometres / KILOMETRES_IN_MILE The first thing to note is that we create and initialize a few variables (or "pseudo-constants" denoted by upper case names) for later use: DEGREES_IN_RADIAN, MEAN_EARTH_RADIUS_KM and KILOMETRES_IN_MILE.
The main functionality is contained in a class called GreatCircle, and in __init__ we create a number of properties with default values. Most pieces of information exist twice: the latitudes, longitudes and angle in both radians and degrees, and the distance in miles and kilometres. There is also a boolean valid flag, a single method called calculate, and a few private functions (denoted with the double underscore or "dunder" prefix).
The greatcircle_calculate function is central to the whole project as it calculates the missing values from the latitudes and longitudes in degrees. Before it does this it checks the degrees are valid with the validate_degrees function, and if so calls calculate_radians, calculate_central_angle and calculate_distance.
Now let's move on to the four private functions called by greatcircle_calculate. The first is validate_degrees which simply checks the latitudes are between -90° and +90°, and that the longitudes are between -180° and +180°.
The calculate_radians function does its stuff by dividing the degrees by DEGREES_IN_RADIAN which was initialised at the top of the file.
The calculate_central_angle function calculates the absolute difference between the longitudes, and then implements this formula:
which we saw above to actually calculate the angle in radians. Finally it calculates the degree equivalent.
The last function in this file is calculate_distance, which as I have mentioned is simply a matter of multiplying the angle in radians by the radius of the earth. We also do the kilometres to miles conversion here.
The class is now finished so let's write a bit of code in main.py to try it out.
import greatcircle
def main():
"""
Creates two lists of dictionaries for start and destination city names and locations.
Iterates list, setting GreatCircle object attributes,
then calculating and printing full details.
"""
print("-------------------------")
print("| codedrome.com |")
print("| Great Circle Distance |")
print("-------------------------\n")
starting_cities = [{"name": "London", "latitude1_degrees": 51.507222, "longitude1_degrees": -0.1275},
{"name": "London", "latitude1_degrees": 51.507222, "longitude1_degrees": -0.1275},
{"name": "London", "latitude1_degrees": 51.507222, "longitude1_degrees": -0.1275},
{"name": "London", "latitude1_degrees": 51.507222, "longitude1_degrees": -0.1275},
{"name": "London", "latitude1_degrees": 51.507222, "longitude1_degrees": -0.1275},
{"name": "London", "latitude1_degrees": 51.507222, "longitude1_degrees": -0.1275},
{"name": "London", "latitude1_degrees": 51.507222, "longitude1_degrees": -0.1275}]
destination_cities = [{"name": "Tokyo", "latitude1_degrees": 35.683333, "longitude1_degrees": 139.683333},
{"name": "New York", "latitude1_degrees": 40.7127, "longitude1_degrees": -74.0059},
{"name": "New Delhi", "latitude1_degrees": 28.613889, "longitude1_degrees": 77.208889},
{"name": "Sydney", "latitude1_degrees": -33.865, "longitude1_degrees": 151.209444},
{"name": "Cape Town", "latitude1_degrees": -33.925278, "longitude1_degrees": 18.423889},
{"name": "Rio de Janeiro", "latitude1_degrees": -22.908333, "longitude1_degrees": -43.196389},
{"name": "Oblivion", "latitude1_degrees": 91, "longitude1_degrees": 360}]
gc = greatcircle.GreatCircle()
for i in range(0, len(starting_cities)):
gc.name1 = starting_cities[i]["name"]
gc.latitude1_degrees = starting_cities[i]["latitude1_degrees"]
gc.longitude1_degrees = starting_cities[i]["longitude1_degrees"]
gc.name2 = destination_cities[i]["name"]
gc.latitude2_degrees = destination_cities[i]["latitude1_degrees"]
gc.longitude2_degrees = destination_cities[i]["longitude1_degrees"]
gc.calculate()
output(gc)
def output(gc):
"""
Prints out all details in a neatly annotated format.
"""
if gc.valid == True:
print("Name 1: " + gc.name1)
print("Latitude %f degrees / %f radians" % (gc.latitude1_degrees, gc.latitude1_radians))
print("Longitude %f degrees / %f radians" % (gc.longitude1_degrees, gc.longitude1_radians))
print("Name 2: " + gc.name2)
print("Latitude %f degrees / %f radians" % (gc.latitude2_degrees, gc.latitude2_radians))
print("Longitude %f degrees / %f radians" % (gc.longitude2_degrees, gc.longitude2_radians))
print("Valid: " + str(gc.valid))
print("Central angle %f radians / %f degrees" % (gc.central_angle_radians, gc.central_angle_degrees))
print("Distance %f kilometers / %f miles" % (gc.distance_kilometres, gc.distance_miles))
else:
print("Latitude and/or longitude for %s and %s are invalid" % (gc.name1, gc.name2))
print("")
main()In the main function I have created two lists of dictionaries representing cities along with their latitudes and longitudes. The last set of values includes the invalid 91° north just to see what happens...!
We then create a GreatCircle object and iterate the lists, setting the object's properties from the current dictionary values before calling calculate. Finally we call the output function which simply prints out the various properties of the object, or a message if the latitude or longitude are invalid.
The code is now finished so we can run it with the following command:
python3 main.py
This is the output - I have just shown the first pair of cities here for brevity.
As well as the names, latitudes and longitudes we supplied we now have the latitudes and longitudes in radians, the central angle in both degrees and radians, and the distance in kilometres and miles.