Sperical (kilometers/miles/radius) distance between two points

gistfile1.py

def earth_radius_at_latitude(lat, unit='km'):
    """Calculates the earth radius at a specific latitude 
    (http://en.wikipedia.org/wiki/Earth_radius)

    Args:
        lat: float representing the latitude of a point
        unit: string (km/miles) representing the measure unit

    Returns:
        the earth radius in the specified unit

    """
    equator_radius = 6378.137
    polar_radius = 6356.7523142
    
    if unit=='miles':
        equator_radius = 3968
        polar_radius = 3947.4
    
    sinlat_pow = math.pow(math.sin(lat), 2)
    
    return equator_radius * math.sqrt(math.pow(polar_radius, 4)/math.pow(equator_radius, 4)*
        sinlat_pow + math.pow(math.cos(lat), 2))/ math.sqrt(1-(1-(math.pow(polar_radius, 2)/math.pow(equator_radius, 2)))*sinlat_pow)
    

def spherical_dist(pos1, pos2, r=3958.75):
    """Calculates the distance between locations, based on each point's longitude and latitude.
    (Harvesine formula)

    Args:
        pos1: array of 1D arrays containing the longitude and latitude of a point.
            example: [9.1, 45.3] or [[9.1, 45.3], [19.3, 20.3]]
        pos2: array of 1D arrays containing the longitude and latitude of a point.
            example: [9.1, 45.3] or [[9.1, 45.3], [19.3, 20.3]]
        r: Multiplicator to trasform the distance from radians to miles or kilometers.
            For example the (mean) Earth radius in miles is 3958.75.
            http://en.wikipedia.org/wiki/Earth_radius

    Returns:
        The distances between the given points in radians(r=1), miles or kilometers.

    """
    pos1 = pos1 * np.pi / 180
    pos2 = pos2 * np.pi / 180
    cos_lat1 = np.cos(pos1[..., 0])
    cos_lat2 = np.cos(pos2[..., 0])
    cos_lat_d = np.cos(pos1[..., 0] - pos2[..., 0])
    cos_lon_d = np.cos(pos1[..., 1] - pos2[..., 1])
    return r * np.arccos(cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))