• Explore
• My scripts
• Calculator

Computer great circle distance between two cities whose coordinates are known

from math import *

# distance using Meeus approximation
def distance( lat1, lon1, lat2, lon2 ) :

    F = (lat1 + lat2)/2.0
    G = (lat1 - lat2)/2.0
    L = (lon1 - lon2)/2.0

    sinG2 = sin(G)**2
    cosG2 = cos(G)**2
    sinF2 = sin(F)**2
    cosF2 = cos(F)**2
    sinL2 = sin(L)**2
    cosL2 = cos(L)**2

    S = sinG2*cosL2 + cosF2*sinL2
    C = cosG2*cosL2 + sinF2*sinL2

    w = atan(sqrt(S/C))
    R = sqrt(S*C)/w

    a = 6378.137            # WGS-84 equatorial radius
    f = 1.0/298.257223563   # WGS-84 ellipsoid flattening factor

    D = 2*w*a
    H1 = (3*R - 1)/(2*C)
    H2 = (3*R + 1)/(2*S)

    dist = D*(1 + f*H1*sinF2*cosG2 - f*H2*cosF2*sinG2)

    return round(100*dist)/100.0


# course using spherical trigonometry
# does not work if one latitude is polar!!!
def course( lat1, lon1, lat2, lon2 ):

    L    = lon2 - lon1
    cosD = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(L)  
    D    = acos(cosD)
    cosC = (sin(lat2) - cosD*sin(lat1))/(sin(D)*cos(lat1))

# numerical error can result in |cosC| slightly > 1.0 
    if cosC > 1.0:
        cosC = 1.0

    if cosC < -1.0:
        cosC = -1.0

    C = 180.0*acos(cosC)/pi

    if sin(L) < 0.0:
        C = 360.0 - C


    return round(100*C)/100.0


# convert to radians
def latitude( sgn, deg, min, sec ):

    lat = 0

    if (0.0 <= sec) and (sec < 60.0):
        lat = lat + sec/3600.0
    else:
        print("check latitude seconds!") 


    if (0.0 <= min) and (min < 60.0):
        lat = lat + min/60.0
    else:
        print("check latitude minutes!")

    
    if (0.0 <= deg) and (deg <= 90.0):
        lat = lat + deg
    else:
        print("check latitude degrees!")


    if (0.0 <= lat) and (lat <= 90.0):
        lat = pi*lat/180.0
    else:
        print("latitude range error!")


    assert (sgn == +1) or (sgn == -1)
    return sgn*lat


# convert to radians
def longitude( sgn, deg, min, sec ):

    lon = 0

    if (0.0 <= sec) and (sec < 60.0):
        lon = lon + sec / 3600.0
    else:
        print("check longitude seconds!")

    
    if (0.0 <= min) and (min < 60.0):
        lon = lon + min / 60.0
    else:
        print("check longitude minutes!")


    if (0.0 <= deg) and (deg <= 180.0):
        lon = lon + deg
    else:
        print("check longitude degrees!")


    if (0.0 <= lon) and (lon <= 180.0):
        lon = pi * lon / 180.0
    else:
        print("longitude range error!")


    assert (sgn == +1) or (sgn == -1)
    return sgn*lon


# Sydney, AS     33°53' S  151°13' E
lat1 = latitude(  -1,  33, 53, 0 )
lon1 = longitude( +1, 151, 13, 0 )

# London, UK,    51°30' N    0°07' W
lat2 = latitude(  +1, 51, 30, 0 )
lon2 = longitude( -1,  0,  7, 0 )

# Boston, MA US  42°21' N   71°04' W
lat3 = latitude(  +1, 42, 21, 0 )
lon3 = longitude( -1, 71,  4, 0 )

# Tokyo, JA,     35°41' N  139°45' E
lat4 = latitude(  +1,  35, 41, 0 )
lon4 = longitude( +1, 139, 45, 0 )

print("Sydney, AS to London, UK")
dist = distance(lat1, lon1, lat2, lon2)
crse = course(lat1, lon1, lat2, lon2)
print("distance = %8.2f  kilometers" %(dist))
print("heading  = %8.2f° true (initially)\n" %(crse))

print("London, UK to Boston, MA US")
dist = distance(lat2, lon2, lat3, lon3)
crse = course(lat2, lon2, lat3, lon3)
print("distance = %8.2f  kilometers" %(dist))
print("heading  = %8.2f° true (initially)\n" %(crse))

print("Boston, MA US to Tokyo, JA")
dist = distance(lat3, lon3, lat4, lon4)
crse = course(lat3, lon3, lat4, lon4)
print("distance = %8.2f  kilometers" %(dist))
print("heading  = %8.2f° true (initially)\n" %(crse))

print("Tokyo, JA to Sydney, AS")
dist = distance(lat4, lon4, lat1, lon1)
crse = course(lat4, lon4, lat1, lon1)
print("distance = %8.2f  kilometers" %(dist))
print("heading  = %8.2f° true (initially)\n" %(crse))