求含经纬度与深度的3D地理坐标间距离的Python实现方法

MATLAB的ecefOffset核心逻辑是将地理坐标（纬度、经度、高度/深度）转换为**地心地固坐标系（ECEF）**的笛卡尔坐标，再计算两点间的三维欧氏距离。以下是对应的Python实现：

实现代码 #

import math

def geodetic_to_ecef(spheroid, lat, lon, h):
    """将地理坐标（纬度、经度、高度/深度）转换为ECEF坐标"""
    # 解析椭球体参数：长半轴a，扁率f
    a = spheroid['a']
    f = spheroid['f']
    b = a * (1 - f)  # 短半轴
    
    # 角度转弧度
    lat_rad = math.radians(lat)
    lon_rad = math.radians(lon)
    
    e_sq = 2 * f - f**2  # 第一偏心率平方
    N = a / math.sqrt(1 - e_sq * math.sin(lat_rad)**2)  # 卯酉圈曲率半径
    
    x = (N + h) * math.cos(lat_rad) * math.cos(lon_rad)
    y = (N + h) * math.cos(lat_rad) * math.sin(lon_rad)
    z = (N * (1 - e_sq) + h) * math.sin(lat_rad)
    
    return x, y, z

def ecef_offset(spheroid, lat1, lon1, h1, lat2, lon2, h2):
    """对应MATLAB ecefOffset，返回两点在ECEF坐标系下的偏移量"""
    x1, y1, z1 = geodetic_to_ecef(spheroid, lat1, lon1, h1)
    x2, y2, z2 = geodetic_to_ecef(spheroid, lat2, lon2, h2)
    
    delta_x = x2 - x1
    delta_y = y2 - y1
    delta_z = z2 - z1
    
    return delta_x, delta_y, delta_z

# 示例用法
if __name__ == "__main__":
    # 使用WGS84椭球体（与MATLAB默认基准一致）
    wgs84_spheroid = {'a': 6378137.0, 'f': 1/298.257223563}
    
    # 注意：深度用负数表示（ECEF高度以椭球面向上为正，深度向下则为负）
    lat1, lon1, depth1 = 39.9042, 116.4074, -100  # 北京附近，深度100米
    lat2, lon2, depth2 = 31.2304, 121.4737, -200  # 上海附近，深度200米
    
    dx, dy, dz = ecef_offset(wgs84_spheroid, lat1, lon1, depth1, lat2, lon2, depth2)
    3d_distance = math.hypot(math.hypot(dx, dy), dz)
    
    print(f"三维空间距离: {3d_distance:.2f} 米")

关键说明 #

1. 椭球体参数：示例用WGS84基准，和MATLAB的referenceEllipsoid('WGS84')完全匹配；若需其他椭球，只需修改spheroid字典的a（长半轴）和f（扁率）。
2. 深度处理：ECEF坐标系中高度以椭球面向上为正，因此深度需传入负数，与MATLAB中h参数的定义一致。
3. 距离计算：通过两次math.hypot计算三维欧氏距离，逻辑和MATLAB的hypot(hypot(deltaX,deltaY),deltaZ)完全等价。

内容的提问来源于stack exchange，提问作者Adnan