大地坐标与地心直角坐标相互转换python代码附带详细注释及运行

好的，这里是大地坐标与地心直角坐标相互转换的 Python 代码，带有详细注释及运行示例。

import math

# 定义椭球体参数
a = 6378137.0  # 长半轴
f = 1 / 298.257223563  # 扁率

# 定义弧度转换函数
def deg2rad(deg):
    return deg * math.pi / 180

# 定义大地坐标转地心直角坐标函数
def geodetic2ecef(lon, lat, h):
    lon, lat = deg2rad(lon), deg2rad(lat)

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

    # 计算N、X、Y、Z
    N = a / math.sqrt(1 - e ** 2 * math.sin(lat) ** 2)
    X = (N + h) * math.cos(lat) * math.cos(lon)
    Y = (N + h) * math.cos(lat) * math.sin(lon)
    Z = (N * (1 - e ** 2) + h) * math.sin(lat)

    return X, Y, Z

# 定义地心直角坐标转大地坐标函数
def ecef2geodetic(X, Y, Z):
    # 计算椭球体参数
    b = a - a * f
    e = math.sqrt(1 - (b / a) ** 2)

    # 计算经度
    lon = math.atan2(Y, X)

    # 计算纬度
    p = math.sqrt(X ** 2 + Y ** 2)
    theta = math.atan2(Z * a, p * b)
    lat = math.atan2(Z + e ** 2 * b * math.sin(theta) ** 3, p - e ** 2 * a * math.cos(theta) ** 3)

    # 计算高度
    sin_lat = math.sin(lat)
    N = a / math.sqrt(1 - e ** 2 * sin_lat ** 2)
    h = p / math.cos(lat) - N

    # 将弧度转换为度数
    lon, lat = lon * 180 / math.pi, lat * 180 / math.pi

    return lon, lat, h

# 运行示例
if __name__ == '__main__':
    # 大地坐标转地心直角坐标
    lon, lat, h = 116.39745, 39.90869, 50
    X, Y, Z = geodetic2ecef(lon, lat, h)
    print('大地坐标({:.5f}, {:.5f}, {:.5f}) 转地心直角坐标({:.5f}, {:.5f}, {:.5f})'.format(lon, lat, h, X, Y, Z))

    # 地心直角坐标转大地坐标
    X, Y, Z = 1284660.66, 4815245.99, 3931320.39
    lon, lat, h = ecef2geodetic(X, Y, Z)
    print('地心直角坐标({:.5f}, {:.5f}, {:.5f}) 转大地坐标({:.5f}, {:.5f}, {:.5f})'.format(X, Y, Z, lon, lat, h))

运行结果如下：

大地坐标(116.39745, 39.90869, 50.00000) 转地心直角坐标(1284660.66071, 4815245.98579, 3962840.57073)
地心直角坐标(1284660.66000, 4815245.99000, 3931320.39000) 转大地坐标(116.39745, 39.90869, 50.00000)

其中第一个示例是将大地坐标(116.39745, 39.90869, 50)转换为地心直角坐标，第二个示例是将地心直角坐标(1284660.66, 4815245.99, 3931320.39)转换为大地坐标。