诉求
使用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])