Python中地理坐标[lon,lat]转保大圆距离笛卡尔[X,Y]坐标方法

不存在可以全局精准保留任意两点间大圆距离的二维平面坐标转换方案。
球面是不可展曲面，根据微分几何的高斯绝妙定理，任何把球面完整摊平到二维平面的操作，必然会在距离、面积、角度三类属性中至少产生一类变形，不可能做到所有点两两之间的平面直线距离，和球面上的大圆最短路径完全相等。

如果你的数据集覆盖范围很小（比如单城市、区域跨度不超过100km），或者只需要保留某一个固定中心点到其余所有点的大圆距离，可以通过成熟的投影算法实现，要么误差小到工程上完全可以忽略，要么特定点对的距离可以做到零偏差。

根据你的数据覆盖范围和精度要求，可以选以下两类落地性强的方案：

• 方位等距投影：选定一个参考点作为投影中心，中心点到平面上任意一点的直线距离，和两点间大圆距离完全相等，中心点到任意点的方位角也完全保真。缺点是两个都不经过中心点的点对，平面距离和真实大圆距离会有偏差，离中心点越远偏差越大，适合所有分析围绕单个核心点展开的场景。
• 局部切平面投影（ENU站心坐标）：取数据集覆盖范围的几何中心作为原点，建立东向、北向的二维切平面坐标系。在区域跨度小于100km的场景下，区域内任意两点的平面欧氏距离和真实大圆距离的相对误差可以低于0.01%，完全满足绝大多数开发场景的精度要求。数据覆盖范围越大误差越高，跨城、跨省、全球尺度的数据集不适用。
• 如果你的需求是全球范围内所有点对的大圆距离都100%保真，不存在对应的二维平面转换方案，只能直接使用三维球面坐标做存储和计算。

直接使用地理空间转换通用库pyproj即可实现上述所有投影转换，不需要自己手写算法逻辑。
首先安装依赖：

pip install pyproj

局部小范围场景（ENU平面坐标）示例代码 #

from pyproj import Transformer
import numpy as np

# 输入坐标格式为[lat, lon]，单位为WGS84坐标系下的度
coords = np.array([
    [39.9042, 116.4074],
    [39.9142, 116.4174],
    [39.8942, 116.3974]
])

# 计算数据集中心点作为投影原点
center_lat = coords[:, 0].mean()
center_lon = coords[:, 1].mean()

# 定义WGS84经纬度转局部平面坐标的转换器
transformer = Transformer.from_proj(
    "EPSG:4326",  # 原始坐标为WGS84经纬度
    f"+proj=aeqd +lat_0={center_lat} +lon_0={center_lon} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs",
    always_xy=True
)

# 转换得到x,y平面坐标，单位为米
# 注意转换器输入顺序为lon, lat，因此对原始坐标做维度调换
x, y = transformer.transform(coords[:, 1], coords[:, 0])
plane_coords = np.column_stack([x, y])
print(plane_coords)

上述代码输出的平面坐标，在小范围场景下两点间的欧氏距离和真实大圆距离的差异可以忽略。

固定中心点的方位等距投影实现 #

如果你只需要保证中心点到所有点的距离完全准确，只需要把上述代码里的center_lat、center_lon改成你选定的固定中心点经纬度即可，不需要取数据集均值，此时中心点到所有其他点的平面欧氏距离和大圆距离完全相等。

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