haversine_numpy.py

"""haversine function adapted to work with numpy arrays

Adapted from https://github.com/mapado/haversine/blob/master/haversine/__init__.py
License: MIT
"""

import math
import numpy as np

AVG_EARTH_RADIUS = 6371 * 1000  # m

def haversine(a, b):
    """ Calculate the great-circle distance between two points on the Earth surface.
    `a` and `b` are numpy arrays with shape (-1, 2). Input array must have the format
    (lon, lat) in degress

    Returns a np.array of distances in meters

    >>> import math
    >>> lyon = (4.8422, 45.7597) 
    >>> paris = (2.3508, 48.8567)
    >>> haversine(np.array([lyon, paris]) , np.array([paris, lyon]))
    array([ 392216.7178066,  392216.7178066])
    """
    radians = math.pi / 180
    lng1, lat1 = a[:, 0] * radians, a[:, 1] * radians
    lng2, lat2 = b[:, 0] * radians, b[:, 1] * radians

    # calculate haversine
    dlat = lat2 - lat1
    dlng = lng2 - lng1
    d = np.sin(dlat * 0.5) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlng * 0.5) ** 2
    h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d))
    return h  # meters