Python实现空间直角坐标转高斯克吕格平面坐标

有时需要将点展开到地图上，需要的是平面坐标，而手中只有三维的空间直角坐标XYZ或者BLH时该怎么办呢？这时候就需要将坐标进行投影转换，本文给出的程序代码是用高斯克吕格投影进行坐标转换

原理示意

源码

源码中还有其他坐标转换的小函数，也可参考

# -*- coding:utf-8 -*-
"""

created on Wed Apr 29 20:30:25 2021
@author: xymeng

"""
import numpy as np
import math as mt
import matplotlib.pyplot as plt
import sklearn.metrics as ms
import csv
import pandas as pd


'''以下为三角函数调用'''
sqrt = mt.sqrt
sin = mt.sin
cos = mt.cos
atan = mt.atan
tan = mt.tan

'''以下为WGS84椭球参数'''
a = 6378137.0         # 参考椭球的长半轴, 单位 m
b = 6356752.31414     # 参考椭球的短半轴, 单位 m
e = sqrt((a**2-b**2)/(a**2)) #椭球第一偏心率
e1 = sqrt((a**2-b**2)/(b**2)) #椭球第二偏心率

'''以下为高斯投影中的参数'''
A1 = 1 + 3/4 * e**2 + 45/64 * e**4 +175/256 * e**6 +11025/16384 * e**8 +43659/65536 * e**10
B1 =     3/4 * e**2 + 15/16 * e**4 +525/512 * e**6 + 2205/2048 * e**8 + 72765/65536 * e**10
C1 =                  15/64 * e**4 +105/256 * e**6 + 2205/4096 * e**8 + 10395/16384 * e**10
D1 =                                35/512 * e**6  + 315/2048 * e**8 + 31185/131072 * e**10
E1 =                                                 315/16384 * e**8  + 3645/65536 * e**10

''''以下为弧度转换角度函数'''
def rad2angle(r):
    """
    该函数可以实现弧度到角度的转换.
    :param r:  弧度
    :return:  a, 对应的角度
    """
    a = r*180.0/mt.pi
    return a

'''以下为角度转换弧度函数'''
def angle2rad(a):
    """
    该函数可以实现角度到弧度的转换.
    :param a:  角度
    :return:  r, 对应的弧度
    """
    r = a*mt.pi/180.0
    return r

'''以下为XYZ转BLH函数'''
def XYZ2BLH(X, Y, Z, a, b):
    """
    该函数实现把某点在参心空间直角坐标系下的坐标（X, Y, Z)转为大地坐标（B, L, H).
    :param X:  X方向坐标，单位 m
    :param Y:  Y方向坐标, 单位 m
    :param Z:  Z方向坐标, 单位 m
    :param a: 地球长半轴，即赤道半径，单位 m
    :param b: 地球短半轴，即大地坐标系原点到两级的距离, 单位 m
    :return:  B, L, H, 大地纬度、经度、海拔高度 (m)
    """

    e = sqrt((a**2-b**2)/(a**2)) #椭球第一偏心率

    if X == 0 and Y > 0:
        L = 90
    elif X == 0 and Y < 0:
        L = -90
    elif X < 0 and Y >= 0:
        L = atan(Y/X)
        L = rad2angle(L)
        L = L+180
    elif X < 0 and Y <= 0:
        L = atan(Y/X)
        L = rad2angle(L)
        L = L-180
    else:
        L = atan(Y / X)
        L = rad2angle(L)

    b0 = atan(Z/sqrt(X**2+Y**2))
    N_temp = a/sqrt((1-e*e*sin(b0)*sin(b0)))
    b1 = atan((Z+N_temp*e*e*sin(b0))/sqrt(X**2+Y**2))

    while abs(b0-b1) > 1e-7:
        b0 = b1
        N_temp = a / sqrt((1 - e * e * sin(b0) * sin(b0)))
        b1 = atan((Z + N_temp * e * e * sin(b0)) / sqrt(X ** 2 + Y ** 2))

    B = b1
    N = a/sqrt((1-e*e*sin(B)*sin(B)))
    H = sqrt(X**2+Y**2)/cos(B)-N
    B = rad2angle(B)
    return B, L, H   # 返回大地纬度B、经度L、海拔高度H (m)

'''以下为高斯克吕格投影6度分带法，经度Lon与投影带号Zone以及带号Zone与投影带中央子午线经度L0的计算'''
def GK(L):
    '''
    高斯克吕格投影6度分带法，经度Lon与投影带号Zone
    以及带号Zone与投影带中央子午线经度L0的计算
    :param L:弧度制
    :return L0:
    '''
    if L>=0:
        '''东半球时'''
        z = int(L/6) + 1
        L0 = (6 * (z) - 3) * mt.pi
    else:
        '''西半球时'''
        z = int(L/6) + 60
        L0 = ((6 * z - 3) - 360)* mt.pi
    return L0

'''以下为将大地坐标转换为高斯平面坐标函数'''
def GKxy(B,L,a,b,A1,B1,C1,D1,E1):
    '''
    将大地坐标转换为高斯平面坐标
    :param B:
    :param L:
    :param a:
    :param b:
    :param A1:
    :param B1:
    :param C1:
    :param D1:
    :param E1:
    :return x,y:
    '''
    X0 = a * (1 - e**2)*(A1 * B - (B1/2)*sin(2*B) + (C1/4)*sin(4*B) - (D1/6)*sin(6*B) + (E1/8)*sin(8*B))
    N = a / sqrt(1 - e*e *sin(B)*sin(B))
    t = tan(B)
    L0 = GK(L)
    l = L - L0
    y2 = e1*e1*cos(B)*cos(B)
    x = X0 + (N/2)*sin(B)*cos(B)*l*l +(N/24)*sin(B)*(cos(B)**3)*(5-t**2 +9*y2 + 4*y2*y2)*(l**4) + (N/720)*sin(B)*(cos(B)**5)*(61 - 58*t*t + (t)**4)*(l**6)
    y = N*cos(B)*l + (N/6)*(cos(B)**3)*(1-t**2 + y2)*(l**3) + (N/120)*(cos(B)**5)*(5-18*(t**2)+(t**4)+14*y2-58*y2*(t**2))*(l**5)
    return x,y

'''以下定义空间直角坐标的存储列表'''
X = []
Y = []
Z = []
'''以下定义大地坐标存储列表'''
B = []
L = []
H = []
'''以下定义高斯平面坐标存储列表'''
x = []
y = []

if __name__ == '__main__':
    path = input('请输入空间直角坐标源文件：') #C://Users//23646//Downloads//数据//IGSNetwork.csv
    with open(path,'r') as f:
        lines = csv.reader(f)
        lines = list(lines)
        for line in lines[1:]:
            X.append(float(line[1]))
            Y.append(float(line[2]))
            Z.append(float(line[3]))

    for i in range(len(X)):
        B0,L0,H0 = XYZ2BLH(X[i],Y[i],Z[i],a,b)
        B.append(B0 * mt.pi/180)
        L.append(L0 * mt.pi/180)
        H.append(H0 * mt.pi/180)

    '''以下是将BLH的弧度数写入到csv文件中'''
    savepath = 'C://Users//23646//Downloads//数据//IGSBLH-rad.csv'
    data = pd.DataFrame({'B':B,'L':L,'H':H})
    data.to_csv(savepath,index=False)

    '''以下为执行大地坐标转高斯平面坐标代码块'''
    for k in range(len(B)):
        x0 ,y0 = GKxy(B[k],L[k],a,b,A1,B1,C1,D1,E1)
        x.append(x0)
        y.append(y0)

    '''以下是将xy的高斯平面坐标写入到csv文件中'''
    savepath1 = 'C://Users//23646//Downloads//数据//IGSxy.csv'
    data = pd.DataFrame({'x': x, 'y': y})
    data.to_csv(savepath1, index=False)