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()