通过Python代码计算两个GPS坐标点间最短路径的实现方法

根据需求场景的不同，分为无约束椭球面最短路径、道路网约束通行最短路径两类实现方案，所有代码基于Python3环境可直接运行。

场景1：无障碍物/路网约束的椭球面最短路径 #

GPS坐标默认采用WGS84椭球坐标系，两点在椭球面上的最短路径为大地线，适用于航空、航海、两点通视测距等无地表通行限制的场景。

• 核心算法：优先采用Karney大地线反算算法，精度可达纳米级，优于传统Vincenty算法；低精度快速估算可采用Haversine球面近似公式，计算误差小于0.5%。
• Haversine近似公式（地球平均半径取$R=6371000\mathrm{m}$）：
  $$a = \sin^2\left(\frac{\Delta\phi}{2}\right) + \cos\phi_1 \cdot \cos\phi_2 \cdot \sin^2\left(\frac{\Delta\lambda}{2}\right)$$
  $$c = 2 \cdot \mathrm{atan2}\left(\sqrt{a}, \sqrt{1-a}\right)$$
  $$d = R \cdot c$$
  其中$\phi$为纬度、$\lambda$为经度，单位为弧度；$d$为两点间球面近似距离。
• Python实现（基于pyproj的高精度大地线计算）：
  首先安装依赖：
  pip install pyproj

from pyproj import Geod
# 初始化WGS84椭球计算实例
geod_calc = Geod(ellps="WGS84")

def geodesic_shortest_path(start_lon: float, start_lat: float, end_lon: float, end_lat: float, interp_num:int=20):
    """
    计算WGS84椭球面上两点间最短大地线路径
    :param start_lon: 起点经度（十进制度）
    :param start_lat: 起点纬度（十进制度）
    :param end_lon: 终点经度（十进制度）
    :param end_lat: 终点纬度（十进制度）
    :param interp_num: 路径插值点数量，数值越大路径越平滑
    :return: 路径总长度（单位：米）、起点出发方位角（单位：度）、路径坐标点列表[(lon, lat)]
    """
    # 反算两点间方位角与总长度
    azi_start, azi_end, total_dist = geod_calc.inv(start_lon, start_lat, end_lon, end_lat)
    # 沿大地线生成插值点
    interp_points = geod_calc.npts(start_lon, start_lat, end_lon, end_lat, interp_num)
    # 拼接起终点形成完整路径
    full_path = [(start_lon, start_lat)] + interp_points + [(end_lon, end_lat)]
    return total_dist, azi_start, full_path

# 调用示例：起点北京天安门，终点北京西站
if __name__ == "__main__":
    dist, azi, path = geodesic_shortest_path(116.397428, 39.90923, 116.32147, 39.89486)
    print(f"大地线最短距离：{dist:.2f}米，起始方位角：{azi:.2f}度")

场景2：道路网约束下的实际通行最短路径 #

适用于驾车、步行、骑行等需要沿可通行道路行驶的导航场景，路径不可穿越建筑、河流、禁行区域。

• 核心逻辑：将路网抽象为带权有向图，节点为道路交叉口，边为路段，权重为路段长度/通行时间，采用Dijkstra、A*、Contraction Hierarchies等图最短路径算法求解。
• Python离线实现（基于OpenStreetMap开源路网数据）：
  首先安装依赖：
  pip install osmnx networkx

import osmnx as ox
import networkx as nx

# 全局配置，开启本地缓存避免重复下载路网
ox.config(use_cache=True, log_console=False)
# 下载目标区域驾驶路网，network_type可替换为walk/bike对应步行/骑行路网
# 示例下载中心点为北京天安门，半径2000米范围的路网
road_graph = ox.graph_from_point((39.90923, 116.397428), dist=2000, network_type="drive")
# 为路网边添加限速、通行时间权重
road_graph = ox.add_edge_speeds(road_graph)
road_graph = ox.add_edge_travel_times(road_graph)

def road_shortest_path(start_lat: float, start_lon: float, end_lat: float, end_lon: float, weight_type:str="length"):
    """
    求解道路网约束下的最短路径
    :param start_lat: 起点纬度（十进制度）
    :param start_lon: 起点经度（十进制度）
    :param end_lat: 终点纬度（十进制度）
    :param end_lon: 终点经度（十进制度）
    :param weight_type: 路径权重类型，length=最短距离，travel_time=最短通行时间
    :return: 路径总权重（单位：米/秒，对应权重类型）、路径坐标点列表[(lon, lat)]
    """
    # 匹配距离起终点最近的路网节点
    origin_node = ox.nearest_nodes(road_graph, start_lon, start_lat)
    dest_node = ox.nearest_nodes(road_graph, end_lon, end_lat)
    # 求解最短路径节点序列
    route_nodes = nx.shortest_path(road_graph, origin_node, dest_node, weight=weight_type)
    total_weight = nx.shortest_path_length(road_graph, origin_node, dest_node, weight=weight_type)
    # 提取路径几何坐标
    route_geom = ox.routing.route_to_gdf(road_graph, route_nodes)
    path_points = []
    for _, row in route_geom.iterrows():
        coords = list(row["geometry"].coords)
        path_points.extend(coords)
    # 路径点去重
    path_points = list(dict.fromkeys(path_points))
    return total_weight, path_points

# 调用示例
if __name__ == "__main__":
    total_len, path = road_shortest_path(39.90923, 116.397428, 39.89486, 116.32147, weight_type="length")
    print(f"道路网最短路径长度：{total_len:.2f}米")

注意事项 #

• 输入GPS坐标时必须严格确认经纬度顺序：pyproj接口采用(经度, 纬度)顺序，osmnx的点匹配接口采用(纬度, 经度)顺序，顺序错误会导致结果偏移数百公里。
• 若使用自有GPS点位数据集批量求解最短路径，可先通过KDTree构建空间索引，预计算点位间的邻接关系与通行权重，构建自定义图结构后调用networkx的最短路径接口即可，无需依赖公开路网数据。
• 长距离跨区域路径计算禁止直接将经纬度投影到平面做欧氏距离计算，会引入不可忽略的距离误差。

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