python不依赖第三方库进行平面坐标和球面坐标的转换

1. 开发需求：前端无法应用python复杂库，将坐标转换用数学函数进行编写

2. 此篇文章借鉴了非常多的资料，代码是借用其他优秀作者的代码，但时间过去太久忘记了代码出处，此文章仅用作自己学习，如有侵权请联系我删除

3. WGS84球面坐标转UTM平面坐标的代码

#WGS84转UTM
import numpy as np
from math import sqrt,cos,sin,tan  #到如所需数学函数库

def WGS_UTM(point: np.array)->(np.array):
    # 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 # 第一偏心率

    # UTM Parameters 参数
    UTM_K0 = 0.9996 # scale factor 比例因子：0°经线投影后和原始经线的比率
    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
    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])

WGS_UTM((120,38))   #结果：array([ 236578.62391666, 4210063.03409706])

4. UTM平面坐标转WGS84坐标代码

#UTM转WGS84 原作者写的太好以至于没有任何修改
import math

'''####################################################
类XYexchangeBL   进行WGS84坐标系下 的XY转为BL  
                 X已东偏500000
                 需要手动输入当地经线 
####################################################'''
x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 扁率

class XYexchangeBL:
    def get_WGS84_af(self):
        '''
         * WGS84
         * 长半轴a=6378137± 2（m）
         * 短半轴b=6356752.3142m
         * 扁率α=1/298.257223563
         * 第一偏心率平方 =0.00669437999013
         * 第二偏心率平方 =0.00673949674223
        '''
        a=6378137.0
        f=1/298.257223563
        return a,f
    
        '''###################
        #函数:Process_Degree(self,dD)
        #输入值 十进制的经纬度(DEG) 
        #输出 度分秒的经纬度(DMS)
        ###################'''

    def XY2LatLon(self,ellipsoid,X, Y, Longitude):
        X=X/0.9996#仅针对UTM投影
        Y=Y/0.9996#针对UTM投影
        N=int(Longitude/6)
        L0=6*(N+1)-3
        #椭圆参数控制
        if(ellipsoid==84):
            a,f=self.get_WGS84_af()

        iPI=0.0174532925199433333333#圆周率/180
        ProjNo=int(X/1000000)
        L0=L0*iPI


        X0=ProjNo*1000000+500000#东偏500000为后续步骤减去做铺垫
        Y0=0
        xval=X-X0
        yval=Y-Y0

        #e2=2*f-f*f#第一偏心率平方
        #e1=(1.0-math.sqrt(1-e2))/(1.0+math.sqrt(1-e2))
        #ee=e2/(1-e2)#第二偏心率平方
        e2=0.00669437999013
        e1=(1.0-math.sqrt(1-e2))/(1.0+math.sqrt(1-e2))
        ee=0.00673949674223


        M=yval
        u=M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256))

        #"\"表示转公式下一行结合在一起
        fai=u+(3*e1/2-27*e1*e1*e1/32)*math.sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*math.sin(4*u)+(151*e1*e1*e1/96)*math.sin(6*u)+(1097*e1*e1*e1*e1/512)*math.sin(8*u)
        C=ee*math.cos(fai)*math.cos(fai)
        T=math.tan(fai)*math.tan(fai)
        NN=a/math.sqrt(1.0-e2*math.sin(fai)*math.sin(fai))
        R=a*(1-e2)/math.sqrt((1-e2*math.sin(fai)*math.sin(fai))*(1-e2*math.sin(fai)*math.sin(fai))*(1-e2*math.sin(fai)*math.sin(fai)))
        D=xval/NN
        #计算经纬度（弧度单位的经纬度）
        longitude1=L0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D*D*D*D*D/120)/math.cos(fai)
        latitude1=fai-(NN*math.tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720)

        #换换为deg
        #longitude=self.Process_Degree(longitude1/iPI)
        #latitude=self.Process_Degree(latitude1/iPI)

        longitude=longitude1/iPI
        latitude=latitude1/iPI
        #return latitude,longitude
        return [latitude,longitude]
  
can=XYexchangeBL()
can.XY2LatLon(84,236578.62391666, 4210063.03409706,120)