如何在地球表面用线段插值以P0为圆心的P1-P2圆弧？

可以通过pyproj、geopy和numpy组合实现球面圆弧的等距插值，核心思路是借助地心坐标系（ECEF）完成球面几何计算，步骤如下：

• 将经纬度转换为地心直角坐标系，简化球面角度与向量运算
• 计算P1、P2相对于P0的方位角，以及两点间的圆心角
• 根据插值段数N拆分圆心角，逐次计算每个插值点的方位角
• 基于P0坐标、固定距离和目标方位角，反推插值点的经纬度
• 将所有点整理为GeoJSON支持的线段格式

以下是完整实现代码：

import geopy.distance
import numpy as np
from pyproj import Transformer

# 初始化经纬度<->地心坐标系（ECEF）转换器
ecef_to_ll = Transformer.from_crs("EPSG:4978", "EPSG:4326", always_xy=True)
ll_to_ecef = Transformer.from_crs("EPSG:4326", "EPSG:4978", always_xy=True)

# 圆心P0坐标
lon0 = 14.265608
lat0 = 50.095867
coords0 = (lat0, lon0)

# 圆弧端点P1坐标
lon1 = 14.475131
lat1 = 50.052842
coords1 = (lat1, lon1)

# 圆弧端点P2坐标
lon2 = 14.416322
lat2 = 49.992767
coords2 = (lat2, lon2)

# 验证P1、P2到P0的距离相等
dist_p0_p1 = geopy.distance.distance(coords1, coords0).nautical
dist_p0_p2 = geopy.distance.distance(coords2, coords0).nautical
print(f"P0到P1距离（海里）: {dist_p0_p1:.2f}")
print(f"P0到P2距离（海里）: {dist_p0_p2:.2f}")

def interpolate_arc(p0, p1, p2, num_segments):
    # 转换所有点到ECEF坐标系
    x0, y0, z0 = ll_to_ecef.transform(p0[1], p0[0])
    x1, y1, z1 = ll_to_ecef.transform(p1[1], p1[0])
    x2, y2, z2 = ll_to_ecef.transform(p2[1], p2[0])
    
    # 计算P1、P2相对于P0的向量
    vec1 = np.array([x1 - x0, y1 - y0, z1 - z0])
    vec2 = np.array([x2 - x0, y2 - y0, z2 - z0])
    
    # 计算两点间的圆心角（弧度）
    dot_product = np.dot(vec1, vec2)
    norm_product = np.linalg.norm(vec1) * np.linalg.norm(vec2)
    central_angle = np.arccos(np.clip(dot_product / norm_product, -1.0, 1.0))
    
    # 确定旋转轴（两向量的叉乘）
    rotation_axis = np.cross(vec1, vec2)
    rotation_axis /= np.linalg.norm(rotation_axis)
    
    # 生成等间隔的旋转角度序列
    angles = np.linspace(0, central_angle, num_segments + 1)
    
    # 初始化插值点列表，先加入起点P1
    interpolated_points = [p1]
    
    # 计算每个中间插值点
    for angle in angles[1:-1]:
        # 用罗德里格斯旋转公式旋转向量，保证到P0的距离不变
        cos_theta = np.cos(angle)
        sin_theta = np.sin(angle)
        rotated_vec = cos_theta * vec1 + sin_theta * np.cross(rotation_axis, vec1) + (1 - cos_theta) * np.dot(rotation_axis, vec1) * rotation_axis
        
        # 转换回经纬度坐标
        x = x0 + rotated_vec[0]
        y = y0 + rotated_vec[1]
        z = z0 + rotated_vec[2]
        lon, lat = ecef_to_ll.transform(x, y, z)
        interpolated_points.append((lat, lon))
    
    # 加入终点P2
    interpolated_points.append(p2)
    
    return interpolated_points

# 插值为5段（生成6个点）
num_segments = 5
arc_points = interpolate_arc(coords0, coords1, coords2, num_segments)

# 输出所有插值点
print("\n插值后的圆弧点（lat, lon）:")
for i, point in enumerate(arc_points):
    print(f"点{i}: {point[0]:.6f}, {point[1]:.6f}")

# 转换为GeoJSON LineString格式
geojson = {
    "type": "Feature",
    "properties": {},
    "geometry": {
        "type": "LineString",
        "coordinates": [[lon, lat] for lat, lon in arc_points]
    }
}
print("\n适配GeoJSON的格式:")
print(geojson)

代码说明：

• 用pyproj处理坐标系转换，避免直接计算球面几何的复杂公式
• 罗德里格斯旋转公式确保每个插值点到P0的距离严格一致
• 最终输出的点序列可直接作为GeoJSON的线段坐标，适配不支持圆弧的地理数据格式

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