地心直角坐标转大地坐标 python代码 用公式进行编写

将地心直角坐标(x, y, z)转换为大地坐标(经度lon, 纬度lat, 高度h)的公式如下：

首先计算出大地坐标系的原点处的地心直角坐标系坐标(x0, y0, z0)，可以采用WGS84椭球体参数：

a = 6378137.0  # WGS84椭球体长半轴
b = 6356752.3142  # WGS84椭球体短半轴

e2 = 1 - (b/a)**2  # 第一偏心率平方
e2_ = (a/b)**2 - 1  # 第二偏心率平方

# 计算大地坐标系原点处的地心直角坐标系坐标
N0 = a / (math.sqrt(1 - e2 * math.sin(0)**2))  # 计算曲率半径
x0 = (N0 + 0) * math.cos(0)  # 原点处地心直角坐标系x坐标
y0 = (N0 + 0) * math.sin(0)  # 原点处地心直角坐标系y坐标
z0 = ((1 - e2) * N0 + 0) * math.sin(0)  # 原点处地心直角坐标系z坐标

接下来，可以通过以下公式计算大地坐标系中的经度和纬度：

p = math.sqrt(x**2 + y**2)  # 计算投影距离
theta = math.atan2(z * a, p * b)  # 计算纬度角theta
lon = math.atan2(y, x)  # 计算经度角lon

最后，可以通过以下公式计算大地高：

N = a / (math.sqrt(1 - e2 * math.sin(theta)**2))  # 计算曲率半径
h = p / math.cos(theta) - N  # 计算高度h

将上述公式整合到Python代码中，可以实现地心直角坐标系到大地坐标系的转换：

import math

def ecef_to_lla(x, y, z):
    # WGS84椭球体参数
    a = 6378137.0
    b = 6356752.3142
    e2 = 1 - (b/a)**2
    e2_ = (a/b)**2 - 1

    # 计算大地坐标系原点处的地心直角坐标系坐标
    N0 = a / (math.sqrt(1 - e2 * math.sin(0)**2))
    x0 = (N0 + 0) * math.cos(0)
    y0 = (N0 + 0) * math.sin(0)
    z0 = ((1 - e2) * N0 + 0) * math.sin(0)

    # 计算大地坐标系中的经度和纬度
    p = math.sqrt(x**2 + y**2)
    theta = math.atan2(z * a, p * b)
    lon = math.atan2(y, x)

    # 计算大地高
    N = a / (math.sqrt(1 - e2 * math.sin(theta)**2))
    h = p / math.cos(theta) - N

    # 转换为角度制
    lon = math.degrees(lon)
    theta = math.degrees(theta)

    # 输出结果
    print("经度：{:.6f}，纬度：{:.6f}，高度：{:.6f}".format(lon, theta, h))

其中，(x, y, z)为输入的地心直角坐标系坐标，输出结果为转换后的大地坐标系中的经度、纬度和高度。