Longitude & Latitude

lon_lat_tool.py

import math


"""
Sample:
lat = 25.0763881
lon = 121.5259891
brng = 45
dist = 600

lat1 = 25.0763881
lon1 = 121.5259891
lat2 = 25.080218035408
lon2 = 121.53019454571


print distVincenty(lat1, lon1, lat2, lon2)

print place_len_to(lat, lon, brng, dist)
"""

def distVincenty(lat1, long1, lat2, long2):
    """
    get distance between 2 points
    """
    a = 6378137
    b = 6356752.314245
    f = 1/298.257223563
    L = math.radians(long2 - long1)
    U1 = math.atan((1 - f) * math.tan(math.radians(lat1)))
    U2 = math.atan((1 - f) * math.tan(math.radians(lat2)))
    sinU1 = math.sin(U1)
    cosU1 = math.cos(U1)
    sinU2 = math.sin(U2)
    cosU2 = math.cos(U2)
    lamb = L
    iterLimit = 100
    while True:
        sinLambda = math.sin(lamb)
        cosLambda = math.cos(lamb)
        sinSigma = math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda))
        if sinSigma == 0 :
            return 0

        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        sigma = math.atan2(sinSigma, cosSigma)
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha * sinAlpha
        if cosSqAlpha == 0:
            cos2SigmaM = 0
        else:
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha

        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lamb
        lamb = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))

        if abs(lamb - lambdaP) > 0.000000000001 and --iterLimit > 0:
            pass
        else:
            break

    if iterLimit == 0 :
        return 0

    uSq = cosSqAlpha * (a * a - b * b) / (b * b)
    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 * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)))

    return b * A * (sigma - deltaSigma)


def place_len_to(lat0, long0, brng, dist):
    """
    get point via brng and dist

    top: brng = 0
    right: brng = 90
    down: brng = 180
    left: brng = 270
    """
    a = 6378137
    b = 6356752.3142
    f = 1 / 298.257223563
    alpha1 = math.radians(brng)
    sinAlpha1 = math.sin(alpha1)
    cosAlpha1 = math.cos(alpha1)
    tanU1 = (1 - f) * math.tan(math.radians(lat0))
    cosU1 = 1 / math.sqrt((1 + tanU1 * tanU1))
    sinU1 = tanU1 * cosU1
    sigma1 = math.atan2(tanU1, cosAlpha1)
    sinAlpha = cosU1 * sinAlpha1
    cosSqAlpha = 1 - sinAlpha * sinAlpha
    uSq = cosSqAlpha * (a * a - b * b) / (b * b)
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
    sigma = dist / (b * A)
    sigmaP = 2 * math.pi
    while abs(sigma - sigmaP) > 0.000000000001 :
        cos2SigmaM = math.cos(2 * sigma1 + sigma)
        sinSigma = math.sin(sigma)
        cosSigma = math.cos(sigma)
        deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)))
        sigmaP = sigma
        sigma = dist / (b * A) + deltaSigma

    tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1
    lat_end = math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, (1 - f) * math.sqrt(sinAlpha * sinAlpha + tmp * tmp))
    lamb = math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1)
    C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
    L = lamb - (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))
    long_end = math.fmod((math.radians(long0) + L + 3 * math.pi), 2 * math.pi) - math.pi

    return (math.degrees(lat_end), math.degrees(long_end))