haversine.py

from enum import Enum
from math import pi
from typing import Union, Tuple
import math


# mean earth radius - https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
_AVG_EARTH_RADIUS_KM = 6371.0088


class Unit(str, Enum):
    """
    Enumeration of supported units.
    The full list can be checked by iterating over the class; e.g.
    the expression `tuple(Unit)`.
    """

    KILOMETERS = 'km'
    METERS = 'm'
    MILES = 'mi'
    NAUTICAL_MILES = 'nmi'
    FEET = 'ft'
    INCHES = 'in'
    RADIANS = 'rad'
    DEGREES = 'deg'


class Direction(float, Enum):
    """
    Enumeration of supported directions.
    The full list can be checked by iterating over the class; e.g.
    the expression `tuple(Direction)`.
    Angles expressed in radians.
    """

    NORTH = 0.0
    NORTHEAST = pi * 0.25
    EAST = pi * 0.5
    SOUTHEAST = pi * 0.75
    SOUTH = pi
    SOUTHWEST = pi * 1.25
    WEST = pi * 1.5
    NORTHWEST = pi * 1.75


# Unit values taken from http://www.unitconversion.org/unit_converter/length.html
_CONVERSIONS = {
    Unit.KILOMETERS:       1.0,
    Unit.METERS:           1000.0,
    Unit.MILES:            0.621371192,
    Unit.NAUTICAL_MILES:   0.539956803,
    Unit.FEET:             3280.839895013,
    Unit.INCHES:           39370.078740158,
    Unit.RADIANS:          1/_AVG_EARTH_RADIUS_KM,
    Unit.DEGREES:          (1/_AVG_EARTH_RADIUS_KM)*(180.0/pi)
}


def get_avg_earth_radius(unit):
    return _AVG_EARTH_RADIUS_KM * _CONVERSIONS[unit]


def _normalize(lat: float, lon: float) -> Tuple[float, float]:
    """
    Normalize point to [-90, 90] latitude and [-180, 180] longitude.
    """
    lat = (lat + 90) % 360 - 90
    if lat > 90:
        lat = 180 - lat
        lon += 180
    lon = (lon + 180) % 360 - 180
    return lat, lon


def _normalize_vector(lat: "numpy.ndarray", lon: "numpy.ndarray") -> Tuple["numpy.ndarray", "numpy.ndarray"]:
    """
    Normalize points to [-90, 90] latitude and [-180, 180] longitude.
    """
    lat = (lat + 90) % 360 - 90
    lon = (lon + 180) % 360 - 180
    wrap = lat > 90
    if numpy.any(wrap):
        lat[wrap] = 180 - lat[wrap]
        lon[wrap] = lon[wrap] % 360 - 180
    return lat, lon


def _ensure_lat_lon(lat: float, lon: float):
    """
    Ensure that the given latitude and longitude have proper values. An exception is raised if they are not.
    """
    if lat < -90 or lat > 90:
        raise ValueError(f"Latitude {lat} is out of range [-90, 90]")
    if lon < -180 or lon > 180:
        raise ValueError(f"Longitude {lon} is out of range [-180, 180]")


def _ensure_lat_lon_vector(lat: "numpy.ndarray", lon: "numpy.ndarray"):
    """
    Ensure that the given latitude and longitude have proper values. An exception is raised if they are not.
    """
    if numpy.abs(lat).max() > 90:
        raise ValueError("Latitude(s) out of range [-90, 90]")
    if numpy.abs(lon).max() > 180:
        raise ValueError("Longitude(s) out of range [-180, 180]")


def _explode_args(f):
    return lambda ops: f(**ops.__dict__)


@_explode_args
def _create_haversine_kernel(*, asin=None, arcsin=None, cos, radians, sin, sqrt, **_):
    asin = asin or arcsin

    def _haversine_kernel(lat1, lng1, lat2, lng2):
        """
        Compute the haversine distance on unit sphere.  Inputs are in degrees,
        either scalars (with ops==math) or arrays (with ops==numpy).
        """
        lat1 = radians(lat1)
        lng1 = radians(lng1)
        lat2 = radians(lat2)
        lng2 = radians(lng2)
        lat = lat2 - lat1
        lng = lng2 - lng1
        d = (sin(lat * 0.5) ** 2
             + cos(lat1) * cos(lat2) * sin(lng * 0.5) ** 2)
        # Note: 2 * atan2(sqrt(d), sqrt(1-d)) is more accurate at
        # large distance (d is close to 1), but also slower.
        return 2 * asin(sqrt(d))
    return _haversine_kernel



@_explode_args
def _create_inverse_haversine_kernel(*, asin=None, arcsin=None, atan2=None, arctan2=None, cos, degrees, radians, sin, sqrt, **_):
    asin = asin or arcsin
    atan2 = atan2 or arctan2

    def _inverse_haversine_kernel(lat, lng, direction, d):
        """
        Compute the inverse haversine on unit sphere.  lat/lng are in degrees,
        direction in radians; all inputs are either scalars (with ops==math) or
        arrays (with ops==numpy).
        """
        lat = radians(lat)
        lng = radians(lng)
        cos_d, sin_d = cos(d), sin(d)
        cos_lat, sin_lat = cos(lat), sin(lat)
        sin_d_cos_lat = sin_d * cos_lat
        return_lat = asin(cos_d * sin_lat + sin_d_cos_lat * cos(direction))
        return_lng = lng + atan2(sin(direction) * sin_d_cos_lat,
                                 cos_d - sin_lat * sin(return_lat))
        return degrees(return_lat), degrees(return_lng)
    return _inverse_haversine_kernel


_haversine_kernel = _create_haversine_kernel(math)
_inverse_haversine_kernel = _create_inverse_haversine_kernel(math)

try:
    import numpy
    has_numpy = True
    _haversine_kernel_vector = _create_haversine_kernel(numpy)
    _inverse_haversine_kernel_vector = _create_inverse_haversine_kernel(numpy)
except ModuleNotFoundError:
    # Import error will be reported in haversine_vector() / inverse_haversine_vector()
    has_numpy = False

try:
    import numba # type: ignore
    if has_numpy:
        _haversine_kernel_vector = numba.vectorize(fastmath=True)(_haversine_kernel_vector)
        # Tuple output is not supported for numba.vectorize. Just jit the numpy version.
        _inverse_haversine_kernel_vector = numba.njit(fastmath=True)(_inverse_haversine_kernel_vector)
    _haversine_kernel = numba.njit(_haversine_kernel)
    _inverse_haversine_kernel = numba.njit(_inverse_haversine_kernel)
except ModuleNotFoundError:
    pass


def haversine(point1, point2, unit=Unit.KILOMETERS, normalize=False, check=True):
    """ Calculate the great-circle distance between two points on the Earth surface.

    Takes two 2-tuples, containing the latitude and longitude of each point in decimal degrees,
    and, optionally, a unit of length.

    :param point1: first point; tuple of (latitude, longitude) in decimal degrees
    :param point2: second point; tuple of (latitude, longitude) in decimal degrees
    :param unit: a member of haversine.Unit, or, equivalently, a string containing the
                 initials of its corresponding unit of measurement (i.e. miles = mi)
                 default 'km' (kilometers).
    :param normalize: if True, normalize the points to [-90, 90] latitude and [-180, 180] longitude.
    :param check: if True, check that points are normalized.

    Example: ``haversine((45.7597, 4.8422), (48.8567, 2.3508), unit=Unit.METERS)``

    Precondition: ``unit`` is a supported unit (supported units are listed in the `Unit` enum)

    :return: the distance between the two points in the requested unit, as a float.

    The default returned unit is kilometers. The default unit can be changed by
    setting the unit parameter to a member of ``haversine.Unit``
    (e.g. ``haversine.Unit.INCHES``), or, equivalently, to a string containing the
    corresponding abbreviation (e.g. 'in'). All available units can be found in the ``Unit`` enum.
    """

    # unpack latitude/longitude
    lat1, lng1 = point1
    lat2, lng2 = point2

    # normalize points or ensure they are proper lat/lon, i.e., in [-90, 90] and [-180, 180]
    if normalize:
        lat1, lng1 = _normalize(lat1, lng1)
        lat2, lng2 = _normalize(lat2, lng2)
    elif check:
        _ensure_lat_lon(lat1, lng1)
        _ensure_lat_lon(lat2, lng2)

    return get_avg_earth_radius(unit) * _haversine_kernel(lat1, lng1, lat2, lng2)


def haversine_vector(array1, array2, unit=Unit.KILOMETERS, comb=False, normalize=False, check=True):
    '''
    The exact same function as "haversine", except that this
    version replaces math functions with numpy functions.
    This may make it slightly slower for computing the haversine
    distance between two points, but is much faster for computing
    the distance between two vectors of points due to vectorization.
    '''
    if not has_numpy:
        raise RuntimeError('Error, unable to import Numpy, '
                           'consider using haversine instead of haversine_vector.')

    # ensure arrays are numpy ndarrays
    if not isinstance(array1, numpy.ndarray):
        array1 = numpy.array(array1)
    if not isinstance(array2, numpy.ndarray):
        array2 = numpy.array(array2)

    # ensure will be able to iterate over rows by adding dimension if needed
    if array1.ndim == 1:
        array1 = numpy.expand_dims(array1, 0)
    if array2.ndim == 1:
        array2 = numpy.expand_dims(array2, 0)

    # Asserts that both arrays have same dimensions if not in combination mode
    if not comb:
        if array1.shape != array2.shape:
            raise IndexError(
                "When not in combination mode, arrays must be of same size. If mode is required, use comb=True as argument.")

    # unpack latitude/longitude
    lat1, lng1 = array1[:, 0], array1[:, 1]
    lat2, lng2 = array2[:, 0], array2[:, 1]

    # normalize points or ensure they are proper lat/lon, i.e., in [-90, 90] and [-180, 180]
    if normalize:
        lat1, lng1 = _normalize_vector(lat1, lng1)
        lat2, lng2 = _normalize_vector(lat2, lng2)
    elif check:
        _ensure_lat_lon_vector(lat1, lng1)
        _ensure_lat_lon_vector(lat2, lng2)

    # If in combination mode, turn coordinates of array1 into column vectors for broadcasting
    if comb:
        lat1 = numpy.expand_dims(lat1, axis=0)
        lng1 = numpy.expand_dims(lng1, axis=0)
        lat2 = numpy.expand_dims(lat2, axis=1)
        lng2 = numpy.expand_dims(lng2, axis=1)

    return get_avg_earth_radius(unit) * _haversine_kernel_vector(lat1, lng1, lat2, lng2)


def inverse_haversine(point, distance, direction: Union[Direction, float], unit=Unit.KILOMETERS):
    lat, lng = point
    r = get_avg_earth_radius(unit)
    return _inverse_haversine_kernel(lat, lng, direction, distance/r)


def inverse_haversine_vector(array, distance, direction, unit=Unit.KILOMETERS):
    if not has_numpy:
        raise RuntimeError('Error, unable to import Numpy, '
                           'consider using inverse_haversine instead of inverse_haversine_vector.')

    # ensure arrays are numpy ndarrays
    array, distance, direction = map(numpy.asarray, (array, distance, direction))

    # ensure will be able to iterate over rows by adding dimension if needed
    if array.ndim == 1:
        array = numpy.expand_dims(array, 0)

    # Asserts that arrays are correctly sized
    if array.ndim != 2 or array.shape[1] != 2 or array.shape[0] != len(distance) or array.shape[0] != len(direction):
        raise IndexError("Arrays must be of same size.")

    # unpack latitude/longitude
    lat, lng = array[:, 0], array[:, 1]

    r = get_avg_earth_radius(unit)
    return _inverse_haversine_kernel_vector(lat, lng, direction, distance/r)