Python经纬度转UTM代码

诉求

使用Python将经纬度转换为UTM坐标，现有代码普遍为matlab实现。
提供一种简单实现，不依赖数学之外的三方库。

实现

import numpy as np
from math import sqrt,cos,sin,tan
def LL2UTM(point: np.array)->(np.array):
    '''
    Convert (long,lat) to UTM coords.

    Param:
        point: np.array(lon,lat)
    Return:
        np.martix([UTMEasting, UTMNorthing])

    NOTE:
    East Longitudes are positive, West longitudes are negative.
    North latitudes are positive, South latitudes are negative.
    '''
    # Math const
    M_PI = 3.14159265358979323846
    DEG_TO_RAD = (M_PI / 180.0)

    # WGS84 Parameters
    WGS84_A = 6378137.0 # major axis
    WGS84_E = 0.0818191908 # first eccentricity

    # UTM Parameters
    UTM_K0 = 0.9996 # scale factor
    UTM_E2 = (WGS84_E * WGS84_E) # e^2

    Long, Lat = point
    # Make sure the longitude is between -180.00 .. 179.9
    LongTemp = (Long + 180) - int((Long + 180) / 360) * 360 - 180
    LatRad = Lat * DEG_TO_RAD
    LongRad = LongTemp * DEG_TO_RAD
    ZoneNumber = int((LongTemp + 180) / 6) + 1

    if (Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0):
        ZoneNumber = 32

    # Special zones for Svalbard
    if (Lat >= 72.0 and Lat < 84.0):
        if (LongTemp >= 0.0 and LongTemp < 9.0):
            ZoneNumber = 31
        elif (LongTemp >= 9.0 and LongTemp < 21.0):
            ZoneNumber = 33
        elif (LongTemp >= 21.0 and LongTemp < 33.0):
            ZoneNumber = 35
        elif (LongTemp >= 33.0 and LongTemp < 42.0):
            ZoneNumber = 37

    # +3 puts origin in middle of zone
    LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3
    LongOriginRad = LongOrigin * DEG_TO_RAD

    # compute the UTM Zone from the latitude and longitude
    #UTMZone = str(ZoneNumber) + UTMLetterDesignator(Lat)

    eccPrimeSquared = (UTM_E2) / (1 - UTM_E2)

    N = WGS84_A / sqrt(1 - UTM_E2 * sin(LatRad) * sin(LatRad))
    T = tan(LatRad) * tan(LatRad)
    C = eccPrimeSquared * cos(LatRad) * cos(LatRad)
    A = cos(LatRad) * (LongRad - LongOriginRad)

    M = WGS84_A * ((1 - UTM_E2 / 4 - 3 * UTM_E2 * UTM_E2 / 64 -
        5 * UTM_E2 * UTM_E2 * UTM_E2 / 256) *
        LatRad - (3 * UTM_E2 / 8 + 3 * UTM_E2 * UTM_E2 / 32 +
        45 * UTM_E2 * UTM_E2 * UTM_E2 / 1024) *
        sin(2 * LatRad) + (15 * UTM_E2 * UTM_E2 / 256 +
        45 * UTM_E2 * UTM_E2 * UTM_E2 / 1024) * sin(4 * LatRad) -
        (35 * UTM_E2 * UTM_E2 * UTM_E2 / 3072) * sin(6 * LatRad))

    UTMEasting = (UTM_K0 * N * (A + (1 - T + C) * A * A * A / 6 +
                (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A *
                A * A * A / 120) + 500000.0)

    UTMNorthing = (UTM_K0 * (M + N * tan(LatRad) *
                (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 +
                (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A *
                A * A * A * A * A / 720)))

    # 10000000 meter offset for southern hemisphere
    if (Lat < 0):
        UTMNorthing += 10000000.0

    return np.array([UTMEasting, UTMNorthing])