空间直角坐标系转大地坐标系的python代码

将空间直角坐标系转换为大地坐标系需要使用大地坐标系参数，包括椭球体参数和基准面参数。以下是Python代码实现：

import math

# 定义基准面参数
a = 6378137  # 长半轴
f = 1 / 298.257223563  # 扁率

# 定义空间直角坐标系坐标
x = 1000000
y = 2000000
z = 3000000

# 计算椭球体参数
b = a * (1 - f)
e2 = (a**2 - b**2) / a**2
e_ = math.sqrt((a**2 - b**2) / b**2)

# 计算大地坐标系参数
p = math.sqrt(x**2 + y**2)
theta = math.atan2(z * a, p * b)

N = a / math.sqrt(1 - e2 * math.sin(theta)**2)
h = p / math.cos(theta) - N

phi = math.atan2(z + e_**2 * N * math.sin(theta), p - e2 * a * math.cos(theta))

# 将弧度转换为度分秒
def rad_to_dms(rad):
    d = math.degrees(rad)
    m = (d - int(d)) * 60
    s = (m - int(m)) * 60
    return (int(d), int(m), s)

# 输出大地坐标系坐标
print("经度：", rad_to_dms(phi))
print("纬度：", rad_to_dms(theta))
print("高度：", h)

其中，`a`为椭球体的长半轴，`f`为扁率。`x`、`y`、`z`为空间直角坐标系的坐标。`b`为椭球体的短半轴，`e2`为第一偏心率的平方，`e_`为第二偏心率。`p`为水平面距离，`theta`为与Z轴的夹角。`N`为卯酉圈曲率半径，`h`为大地高。`phi`为经度，`theta`为纬度，均以弧度表示，需要用`rad_to_dms`函数将其转换为度分秒形式输出。