def VincentyDistance(lon1_in,lat1_in,lon2_in,lat2_in,ELLIPSOID='GRS-80'):
    """
    Calculate the geodesic distance between two points using the formula
    devised by Thaddeus Vincenty, with an accurate ellipsoidal model of the
    earth.
   
    The class attribute `ELLIPSOID` indicates which ellipsoidal model of the
    earth to use. If it is a string, it is looked up in the `ELLIPSOIDS`
    dictionary to obtain the major and minor semiaxes and the flattening.
    Otherwise, it should be a tuple with those values. The most globally
    accurate model is WGS-84. See the comments above the `ELLIPSOIDS`
    dictionary for more information.
   
    """
    lng1 = np.deg2rad(lon1_in)
    lat1 = np.deg2rad(lat1_in)
    lng2 = np.deg2rad(lon2_in)
    lat2 = np.deg2rad(lat2_in)

    major, minor, f = ELLIPSOIDS[ELLIPSOID]

    delta_lng = lng2 - lng1

    reduced_lat1 = np.arctan((1 - f) * np.tan(lat1))
    reduced_lat2 = np.arctan((1 - f) * np.tan(lat2))

    sin_reduced1, cos_reduced1 = np.sin(reduced_lat1), np.cos(reduced_lat1)
    sin_reduced2, cos_reduced2 = np.sin(reduced_lat2), np.cos(reduced_lat2)

    lambda_lng = delta_lng
    lambda_prime = 2 * np.pi

    iter_limit = 20

    while abs(lambda_lng - lambda_prime) > 10e-12 and iter_limit > 0:
        sin_lambda_lng, cos_lambda_lng = np.sin(lambda_lng), np.cos(lambda_lng)

        sin_sigma = np.sqrt( (cos_reduced2 * sin_lambda_lng) ** 2 +
                             (cos_reduced1 * sin_reduced2 -
                              sin_reduced1 * cos_reduced2 * cos_lambda_lng) ** 2
                            )

        if sin_sigma == 0:
            return 0 # Coincident points

        cos_sigma = (
            sin_reduced1 * sin_reduced2 +
            cos_reduced1 * cos_reduced2 * cos_lambda_lng
        )

        sigma = np.arctan2(sin_sigma, cos_sigma)

        sin_alpha = (
            cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma
        )
        cos_sq_alpha = 1 - sin_alpha ** 2

        if cos_sq_alpha != 0:
            cos2_sigma_m = cos_sigma - 2 * (
                sin_reduced1 * sin_reduced2 / cos_sq_alpha
            )
        else:
            cos2_sigma_m = 0.0 # Equatorial line

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

        lambda_prime = lambda_lng
        lambda_lng = (
            delta_lng + (1 - C) * f * sin_alpha * (
                sigma + C * sin_sigma * (
                    cos2_sigma_m + C * cos_sigma * (
                        -1 + 2 * cos2_sigma_m ** 2
                    )
                )
            )
        )
        iter_limit -= 1

    if iter_limit == 0:
        raise ValueError("Vincenty formula failed to converge!")

    u_sq = cos_sq_alpha * (major ** 2 - minor ** 2) / minor ** 2

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

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

    delta_sigma = (
        B * sin_sigma * (
            cos2_sigma_m + B / 4. * (
                cos_sigma * (
                    -1 + 2 * cos2_sigma_m ** 2
                ) - B / 6. * cos2_sigma_m * (
                    -3 + 4 * sin_sigma ** 2
                ) * (
                    -3 + 4 * cos2_sigma_m ** 2
                )
            )
        )
    )

    s = minor * A * (sigma - delta_sigma)
    return s
def VincentyDistance ( lon1_in , lat1_in , lon2_in , lat2_in , ELLIPSOID = 'GRS-80' ) : lng1 = np . deg2rad ( lon1_in ) lat1 = np . deg2rad ( lat1_in ) lng2 = np . deg2rad ( lon2_in ) lat2 = np . deg2rad ( lat2_in ) major , minor , f = ELLIPSOIDS [ ELLIPSOID ] delta_lng = lng2 - lng1 reduced_lat1 = np . arctan ( ( 1 - f ) * np . tan ( lat1 ) ) reduced_lat2 = np . arctan ( ( 1 - f ) * np . tan ( lat2 ) ) sin_reduced1 , cos_reduced1 = np . sin ( reduced_lat1 ) , np . cos ( reduced_lat1 ) sin_reduced2 , cos_reduced2 = np . sin ( reduced_lat2 ) , np . cos ( reduced_lat2 ) lambda_lng = delta_lng lambda_prime = 2 * np . pi iter_limit = 20 while abs ( lambda_lng - lambda_prime ) > 10e-12 and iter_limit > 0 : sin_lambda_lng , cos_lambda_lng = np . sin ( lambda_lng ) , np . cos ( lambda_lng ) sin_sigma = np . sqrt ( ( cos_reduced2 * sin_lambda_lng ) ** 2 + ( cos_reduced1 * sin_reduced2 - sin_reduced1 * cos_reduced2 * cos_lambda_lng ) ** 2 ) if sin_sigma == 0 : return 0 cos_sigma = ( sin_reduced1 * sin_reduced2 + cos_reduced1 * cos_reduced2 * cos_lambda_lng ) sigma = np . arctan2 ( sin_sigma , cos_sigma ) sin_alpha = ( cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma ) cos_sq_alpha = 1 - sin_alpha ** 2 if cos_sq_alpha != 0 : cos2_sigma_m = cos_sigma - 2 * ( sin_reduced1 * sin_reduced2 / cos_sq_alpha ) else : cos2_sigma_m = 0.0 C = f / 16. * cos_sq_alpha * ( 4 + f * ( 4 - 3 * cos_sq_alpha ) ) lambda_prime = lambda_lng lambda_lng = ( delta_lng + ( 1 - C ) * f * sin_alpha * ( sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * ( - 1 + 2 * cos2_sigma_m ** 2 ) ) ) ) iter_limit -= 1 if iter_limit == 0 : raise ValueError ( "Vincenty formula failed to converge!" ) u_sq = cos_sq_alpha * ( major ** 2 - minor ** 2 ) / minor ** 2 A = 1 + u_sq / 16384. * ( 4096 + u_sq * ( - 768 + u_sq * ( 320 - 175 * u_sq ) ) ) B = u_sq / 1024. * ( 256 + u_sq * ( - 128 + u_sq * ( 74 - 47 * u_sq ) ) ) delta_sigma = ( B * sin_sigma * ( cos2_sigma_m + B / 4. * ( cos_sigma * ( - 1 + 2 * cos2_sigma_m ** 2 ) - B / 6. * cos2_sigma_m * ( - 3 + 4 * sin_sigma ** 2 ) * ( - 3 + 4 * cos2_sigma_m ** 2 ) ) ) ) s = minor * A * ( sigma - delta_sigma ) return s
Calculate the geodesic distance between two points using the formula
    devised by Thaddeus Vincenty, with an accurate ellipsoidal model of the
    earth.
   
    The class attribute `ELLIPSOID` indicates which ellipsoidal model of the
    earth to use. If it is a string, it is looked up in the `ELLIPSOIDS`
    dictionary to obtain the major and minor semiaxes and the flattening.
    Otherwise, it should be a tuple with those values. The most globally
    accurate model is WGS-84. See the comments above the `ELLIPSOIDS`
    dictionary for more information.