Calculating distances is more challenging than compass directions. There are three methods. One is the geometric method which is accurate enough but needs to be recalculated every 50 miles for real accuracy. Another is the Haversine method which is accurate to within about 10% and which was used for well over a hundred years. The best, however is the Vincenty method which is accurate to within 0.5mm.

Here's code comparing all 3 methods...

from vincenty import vincenty

import math

# Python 3 program for the

# haversine formula

def haversine(lat1, lon1, lat2, lon2):

# distance between latitudes

# and longitudes

dLat = (lat2 - lat1) * math.pi / 180.0

dLon = (lon2 - lon1) * math.pi / 180.0

# convert to radians

lat1 = (lat1) * math.pi / 180.0

lat2 = (lat2) * math.pi / 180.0

# apply formulae

a = (pow(math.sin(dLat / 2), 2) +

pow(math.sin(dLon / 2), 2) *

math.cos(lat1) * math.cos(lat2));

rad = 6371

c = 2 * math.asin(math.sqrt(a))

return rad * c

# Driver code

if __name__ == "__main__":

lat1 = 51.5081

lon1 = 0.0759

lat2 = 48.8738

lon2 = -2.2950

print("Haversine distance in miles ",haversine(lat1, lon1,lat2, lon2)/1.6)

# This code is contributed

# by ChitraNayal

LocationOneV = lon1

LocationTwoV = lon2

LocationOneH = lat1

LocationTwoH = lat2

#My own version of distance created using the flat earth model

#Distance calculation

DistanceV = (LocationOneV - LocationTwoV) *54.6

DistanceH = (LocationOneH - LocationTwoH) *68.7

DistanceCalc = (DistanceV*DistanceV)+(DistanceH*DistanceH)

ActualDistance = math.sqrt(DistanceCalc)

print("Distance in Miles using flat earth calculation ", ActualDistance)

TowerOfLondon = (51.5081, 0.0759)

ArcDeTriomph = (48.873756, 2.294946)

print("Vincenty miles ", vincenty(TowerOfLondon, ArcDeTriomph , miles=True))

This, of course, required the Vincenty library which is here:

import math

# WGS 84

a = 6378137 # meters

f = 1 / 298.257223563

b = 6356752.314245 # meters; b = (1 - f)a

MILES_PER_KILOMETER = 0.621371

MAX_ITERATIONS = 200

CONVERGENCE_THRESHOLD = 1e-12 # .000,000,000,001

def vincenty_inverse(point1, point2, miles=False):

"""

Vincenty's formula (inverse method) to calculate the distance (in

kilometers or miles) between two points on the surface of a spheroid

Doctests:

>>> vincenty((0.0, 0.0), (0.0, 0.0)) # coincident points

0.0

>>> vincenty((0.0, 0.0), (0.0, 1.0))

111.319491

>>> vincenty((0.0, 0.0), (1.0, 0.0))

110.574389

>>> vincenty((0.0, 0.0), (0.5, 179.5)) # slow convergence

19936.288579

>>> vincenty((0.0, 0.0), (0.5, 179.7)) # failure to converge

>>> boston = (42.3541165, -71.0693514)

>>> newyork = (40.7791472, -73.9680804)

>>> vincenty(boston, newyork)

298.396057

>>> vincenty(boston, newyork, miles=True)

185.414657

"""

# short-circuit coincident points

if point1[0] == point2[0] and point1[1] == point2[1]:

return 0.0

U1 = math.atan((1 - f) * math.tan(math.radians(point1[0])))

U2 = math.atan((1 - f) * math.tan(math.radians(point2[0])))

L = math.radians(point2[1] - point1[1])

Lambda = L

sinU1 = math.sin(U1)

cosU1 = math.cos(U1)

sinU2 = math.sin(U2)

cosU2 = math.cos(U2)

for iteration in range(MAX_ITERATIONS):

sinLambda = math.sin(Lambda)

cosLambda = math.cos(Lambda)

sinSigma = math.sqrt((cosU2 * sinLambda) ** 2 +

(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ** 2)

if sinSigma == 0:

return 0.0 # coincident points

cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda

sigma = math.atan2(sinSigma, cosSigma)

sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma

cosSqAlpha = 1 - sinAlpha ** 2

try:

cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha

except ZeroDivisionError:

cos2SigmaM = 0

C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))

LambdaPrev = Lambda

Lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma *

(cos2SigmaM + C * cosSigma *

(-1 + 2 * cos2SigmaM ** 2)))

if abs(Lambda - LambdaPrev) < CONVERGENCE_THRESHOLD:

break # successful convergence

else:

return None # failure to converge

uSq = cosSqAlpha * (a ** 2 - b ** 2) / (b ** 2)

A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))

B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))

deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma *

(-1 + 2 * cos2SigmaM ** 2) - B / 6 * cos2SigmaM *

(-3 + 4 * sinSigma ** 2) * (-3 + 4 * cos2SigmaM ** 2)))

s = b * A * (sigma - deltaSigma)

s /= 1000 # meters to kilometers

if miles:

s *= MILES_PER_KILOMETER # kilometers to miles

return round(s, 6)

vincenty = vincenty_inverse

if __name__ == '__main__':

import doctest

doctest.testmod()