Geo Calculations for Flight Navigation

gistfile1.py

from math import radians, cos, sin, asin, sqrt, atan2, degrees

def _c(old):
    direction = {'N':1, 'S':-1, 'E': 1, 'W':-1}
    new = old
    new = new.split()
    new_dir = new.pop(0)
    new.extend([0,0,0])
    x = (float(new[0])+float(new[1])/60.0+float(new[2])/3600.0) * direction[new_dir]
    return x

def bearing(lon1, lat1, lon2, lat2, variation = None):
    if variation:
        variation = variation.split()
        if variation[1].upper() == 'E':
            variation = float(variation[0]) * -1
        else:
            variation = float(variation[0])
    else:
        variation = 0.0
        
    #lon1, lat1, lon2, lat2 = map(radians, [_c(lon1), _c(lat1), _c(lon2), _c(lat2)])
    lon1, lat1, lon2, lat2 = (_c(lon1), _c(lat1), _c(lon2), _c(lat2))

    diffLong = radians(lon2 - lon1)
     
    x = sin(diffLong) * cos(lat2)
    y = cos(lat1) * sin(lat2) - (sin(lat1) * cos(lat2) * cos(diffLong))
     
    initial_bearing = atan2(x, y)
     
    initial_bearing = degrees(initial_bearing)
    compass_bearing = (initial_bearing + 360) % 360
     
    return compass_bearing + variation

def distance(lon1, lat1, lon2, lat2):
    lon1, lat1, lon2, lat2 = map(radians, [_c(lon1), _c(lat1), _c(lon2), _c(lat2)])
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    
    r = 3443.89849 # Radius of earth in NM.
    return c * r

def eet(gs, dist):
    return 60.0 / gs * dist

dist = distance(lon1=u'N 47 27 59.60',lat1=u'E 7 39 55.6',lon2=u'N 47 30 32',lat2=u'E 7 57 0')
bearing = bearing(lon1=u'N 47 27 59.60',lat1=u'E 7 39 55.6',lon2=u'N 47 30 32',lat2=u'E 7 57 0', variation='2 W')
eet_min = eet(100, dist)


print "MC = %.1f / Dist = %.1f / EET = %.1f min" % (bearing, dist, eet_min)