stereographic projections; convert longitude and latitude to meters;

stereographic.py

from typing import (Union,
                    Optional,
                    Tuple)

import numpy as np


def projection(latitudes: np.ndarray,
               longitudes: np.ndarray,
               *,
               latitude_origin: Optional[float] = None,
               longitude_origin: Optional[float] = None,
               earth_radius: float = 6378137,
               earth_eccentricity: float = 0.0818191908426
               ) -> Tuple[np.ndarray, np.ndarray]:
    """
    Transforms latitude/longitude data to map coordinates
    for a polar stereographic system.
    Equations from:
    Map Projections - A Working manual - by J.P. Snyder. 1987
    http://kartoweb.itc.nl/geometrics/Publications/Map%20Projections%20-%20A%20Working%20manual%20-%20by%20J.P.%20Snyder.pdf
    See the section on Stereographic Projection,
    Formulas for the ellipsoid.
    Default Earth radius and eccentricity values are from WGS84.
    :param latitudes: in radians, negative numbers for West
    :param longitudes: in radians, negative numbers for South
    :param latitude_origin: in radians, coordinates origin
    :param longitude_origin: in radians, coordinates origin
    :param earth_radius: radius of Earth in meters
    :param earth_eccentricity: Earth's misshapenness
    :return: X, Y coordinates in meters
    """
    e_sin_phi = get_e_sin_phi(latitude=latitudes,
                              earth_eccentricity=earth_eccentricity)
    chi = get_chi(latitude=latitudes,
                  e_sin_phi=e_sin_phi,
                  earth_eccentricity=earth_eccentricity)

    sin_chi = np.sin(chi)
    cos_chi = np.cos(chi)

    if longitude_origin is None:
        longitude_origin = longitudes[0]

    if latitude_origin is None:
        latitude_origin = latitudes[0]
        sin_chi_1 = sin_chi[0]
        cos_chi_1 = cos_chi[0]
        e_sin_phi_1 = e_sin_phi[0]
    else:
        e_sin_phi_1 = get_e_sin_phi(latitude=latitude_origin,
                                    earth_eccentricity=earth_eccentricity)
        chi_1 = get_chi(latitude=latitude_origin,
                        e_sin_phi=e_sin_phi_1,
                        earth_eccentricity=earth_eccentricity)
        sin_chi_1 = np.sin(chi_1)
        cos_chi_1 = np.cos(chi_1)

    centered_longitudes = longitudes - longitude_origin
    cos_lambda = np.cos(centered_longitudes)

    mu_1 = (np.cos(latitude_origin) / np.sqrt(1 - e_sin_phi_1 ** 2))

    alpha = (2 * earth_radius * mu_1
             / (cos_chi_1 * (1 + sin_chi_1 * sin_chi
                             + cos_chi_1 * cos_chi * cos_lambda)))

    x_values = alpha * cos_chi * np.sin(centered_longitudes)
    y_values = alpha * (cos_chi_1 * sin_chi - sin_chi_1 * cos_chi * cos_lambda)
    return x_values, y_values

def get_e_sin_phi(*,
                  latitude: Union[np.ndarray, float],
                  earth_eccentricity: float) -> Union[np.ndarray, float]:
    return earth_eccentricity * np.sin(latitude)


def get_chi(*,
            latitude: Union[np.ndarray, float],
            e_sin_phi: Union[np.ndarray, float],
            earth_eccentricity: float) -> Union[np.ndarray, float]:
    return -np.pi / 2 + 2 * np.arctan(
        np.tan(np.pi / 4 + latitude / 2)
        * ((1 - e_sin_phi) / (1 + e_sin_phi)) ** (earth_eccentricity / 2))