cheap-ruler: Fast approximations for common geodesic measurements 🌐 (python version)

cheap-ruler.py

import numpy as np

def cheap_ruler_k(latitude: float) -> np.ndarray:
    '''
    in meters
    '''
    # based on https://github.com/mapbox/cheap-ruler-cpp
    RE = 6378.137
    FE = 1.0 / 298.257223563
    E2 = FE * (2 - FE)
    MUL = np.deg2rad(RE * 1000.)
    coslat = np.cos(np.deg2rad(latitude))
    w2 = 1.0 / (1.0 - E2 * (1.0 - coslat * coslat))
    w = np.sqrt(w2)
    k = np.array([MUL * w * coslat, MUL * w * w2 * (1.0 - E2), 1.0])
    k.flags.writeable = False
    return k


def lla2enu(llas: np.ndarray, *,
            anchor_lla: Optional[np.ndarray] = None) -> np.ndarray:
    '''
    WGS84 to ENU
    '''
    llas = np.asarray(llas)
    if llas.shape[1] == 2:
        llas = np.c_[llas, np.zeros(len(llas))]
    if anchor_lla is None:
        anchor_lla = llas[0]
    if len(anchor_lla) == 2:
        anchor_lla = [*anchor_lla, 0.0]
    k = cheap_ruler_k(anchor_lla[1])
    return (llas - anchor_lla) * k


def enu2lla(enus: np.ndarray, *, anchor_lla: np.ndarray) -> np.ndarray:
    '''
    ENU to WGS84
    '''
    enus = np.asarray(enus)
    if enus.shape[1] == 2:
        enus = np.c_[enus, np.zeros(len(enus))]
    k = cheap_ruler_k(anchor_lla[1])
    return anchor_lla + enus / k