transverse_mercator.py

"""
Conversion from WGS-84 to Transverse Mercator
WGS-84 is World Geodetic System data from 1984, utilized by GPS standard.
Transverse mercator is a general Mercator projection (projection of the oblate spheroid of earth onto a cylinder) that is used in most mapping systems. 
These short codes allow user to specify center latitude and longitude.

Contains a simpler projection assuming earth is a sphere
>>> tmercator_meters_to_gps(coords, center)
>>> tmercator_gps_to_meters(coords, center)

and a more precise approximiation leveraging PROJ which ships with current linux systems
proj is a complete Geodetic package https://github.com/OSGeo/PROJ/src/projections/tmerc.cpp
with python bindings https://pypi.org/project/pyproj/

$ pip install pyproj
>>> tmerc(lon, lat, lon_0, lat_0)

Transverse Mercator is sphere to cylinder projection explained in
https://mathworld.wolfram.com/MercatorProjection.html
https://en.wikipedia.org/wiki/Transverse_Mercator_projection

The more precise, yet approximate Evenden Snyder approximation 
https://www.linz.govt.nz/guidance/geodetic-system/understanding-coordinate-conversions/projection-conversions/transverse-mercator-transformation-formulae
"""

from typing import Optional
import numpy as np
from pyproj import Proj

def tmerc(lon, lat, lon_0=None, lat_0=None, inverse=False) -> tuple:
    """ converts from geographic to transverse mercator or back
    centered in lon_0, lat_0
    Args:
        lon     (float, ndarray) in degrees (or meters if inverse)
        lat     (float, ndarray) in degrees (or meters if inverse)
        lon_0   (float) in degrees, center of projection
        lat_0   (float) in degrees
        inverse (bool [False]) if True, meters -> degrees
    """
    lon_0 = lon.min() if lon_0 is None else lon_0
    lat_0 = lat.min() if lat_0 is None else lat_0
    _tmerc = Proj(proj='tmerc', lat_0=lat_0, lon_0=lon_0, ellps='WGS84')
    return _tmerc(lon, lat, inverse=inverse)

###
#
# simple approximation to Transverse Mecator projection assuming earth is a sphere
# https://mathworld.wolfram.com/MercatorProjection.html
# may be ok for some uses. error is ~1.5cm every 10 meters from center
#
def tmercator_gps_meters(coords: np.ndarray,
                         center: Optional [np.ndarray] = None,
                         radius: float = 6_378_137) -> np.ndarray:
    """ degrees to meters using simplified equatorial radius

    Args
        coords  ndarray (N,2) (longitudes, latitudes)
        center  ndarray (2) long, lat, if None, centers on coords
        radius  float default: simplified, equatorial radius
    """
    coords = coords*np.pi/180
    if center is None:
        center = coords.mean(axis=-2)
    center = center*np.pi/180
    long = coords[..., 0]
    lat = coords[..., 1] - center[..., 0]
    # B = cosφ *sin(λ-λ0), x = ln((1+B)/(1-B))/2
    x = np.tanh(np.cos(lat) * np.sin(long))
    y = np.arctan(np.tan(lat)/np.cos(long)) - center[..., 1]
    return radius * np.stack((x, y), axis=-1)

def tmercator_meters_to_gps(pos: np.ndarray,
                            center: np.ndarray,
                            radius: float = 6_378_137) -> np.ndarray:
    """ meters to degrees using simplified equatorial radius

    Args
        pos     ndarray (N,2) (x, y) East, North
        center  ndarray (2) long, lat
        radius  float default: simplified, equatorial radius
    """
    x = pos[...,0] / radius
    y = pos[...,1] / radius
    center = center*np.pi/180
    long = center[..., 0] + np.arctan(np.sinh(x)/np.cos(y + center[..., 1]))
    lat = np.arcsin(np.sin(y + center[..., 1])/np.cosh(x))
    return np.stack((long, lat), axis=-1) * 180 /np.pi