Python中经纬度距离计算的实现方式_Python

python中计算经纬度距离可以使用geopy（vincenty、karney）、haversine、pyproj（椭球模型投影计算）、平面近似法（小范围快速估算）等。而针对中国加密坐标系（如gcj-02/bd-09）也可以使用coord-convert工具转换为wgs84。

以下是对应方法的示例。

一、基本方法

1. 使用geopy库（推荐）

geopy 提供了多种距离计算方式（包括haversine、vincenty、karney），支持wgs84椭球模型。

安装

pip install geopy

示例代码

from geopy.distance import geodesic, great_circle

# 定义两点（纬度, 经度）
point_a = (39.9042, 116.4074)  # 北京
point_b = (31.2304, 121.4737)  # 上海

# 方法1：
# geopy 1.x 版本 vincenty（默认，高精度）
# geopy 2.0 版本 karney
distance_karney = geodesic(point_a, point_b).km
print(f"karney距离: {distance_karney:.2f} km")

# 方法2：大圆距离（haversine，球面近似）
distance_great_circle = great_circle(point_a, point_b).km
print(f"大圆距离: {distance_great_circle:.2f} km")

输出：

karney距离: 1065.85 km
大圆距离: 1067.31 km

说明：
geodesic 2.0版本使用的是karney算法，1.x 版本使用的是vincenty算法。
great_circle 使用haversine公式，假设地球为完美球体。实际地球是椭球体，其误差范围（通常为 0.5% 左右）。

2. 手动实现 haversine 公式

适合无需外部库的场景。

示例代码

import math


def haversine(lat1, lon1, lat2, lon2):
    r = 6371  # 地球半径（km）
    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    delta_phi = math.radians(lat2 - lat1)
    delta_lambda = math.radians(lon2 - lon1)

    a = math.sin(delta_phi / 2) ** 2 + \
        math.cos(phi1) * math.cos(phi2) * math.sin(delta_lambda / 2) ** 2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    return r * c

distance = haversine(39.9042, 116.4074, 31.2304, 121.4737)
print(f"haversine距离: {distance:.2f} km")

输出：

haversine距离: 1067.31 km

注意：与geopy的great_circle结果一致。

3. 使用pyproj进行投影坐标系计算

适合需要平面坐标的高精度场景（如utm投影）。

安装

pip install pyproj

示例代码

from pyproj import geod

# 使用wgs84椭球
geod = geod(ellps="wgs84")

# 计算距离和方位角
# 参数顺序：起点经度、起点纬度、终点经度、终点纬度
_, _, distance = geod.inv(116.4074, 39.9042, 121.4737, 31.2304)
print(f"pyproj椭球距离: {distance / 1000:.2f} km")

输出：

pyproj椭球距离: 1065.85 km

说明：

geod.inv 返回前两个值是方位角，第三个是距离（米）。
结果与geopy的vincenty一致。

4. 平面近似法（小范围适用）

适用于城市内短距离快速估算。

import math

def flat_approximation(lat1, lon1, lat2, lon2):
    r = 6371  # 地球半径（km）
    avg_lat = math.radians((lat1 + lat2) / 2)
    dx = math.radians(lon2 - lon1) * math.cos(avg_lat)
    dy = math.radians(lat2 - lat1)
    return r * math.sqrt(dx**2 + dy**2)

# 示例（北京天安门到北京西站）
distance = flat_approximation(39.9042, 116.4074, 39.8946, 116.3224)
print(f"平面近似距离: {distance:.2f} km")

输出：

平面近似距离: 7.33 km

注意：适用于距离 < 10 km，实际误差还取决于纬度，高纬度地区误差更大。

二、不同坐标系的影响

常见坐标系

许多国家会对地理信息进行加密或偏移，如中国的gcj-02或bd-09，计算距离前需确保所有坐标点在同一坐标系下，否则结果会有偏差。

可以使用库（如coord-convert）将gcj-02或bd-09坐标转换为wgs84，但不同坐标系之间的转换是不可逆的，多次转换可能会导致精度损失。

wgs84（国际标准gps坐标系）：未加密的地心坐标系，高精度。
gcj-02（火星坐标系）：基于wgs84通过非线性算法加入随机偏移。
bd-09（百度坐标系）：在gcj-02基础上二次加密的坐标系。

coord-convert坐标系转换

安装

pip install coord-convert

示例代码

from coord_convert import transform

# gcj-02 转 wgs84
gcj_lon, gcj_lat = 116.404, 39.915  # 北京天安门（gcj-02）
wgs_lon, wgs_lat = transform.gcj2wgs(gcj_lon, gcj_lat)
print(f"wgs84坐标: {wgs_lon:.6f}, {wgs_lat:.6f}")
# wgs84坐标: 116.397756, 39.913596

# bd-09 转 wgs84
bd_lon, bd_lat = 116.410, 39.921  # 北京天安门（bd-09）
wgs_lon, wgs_lat = transform.bd2wgs(bd_lon, bd_lat)
print(f"wgs84坐标: {wgs_lon:.6f}, {wgs_lat:.6f}")
# wgs84坐标: 116.397387, 39.913258

坐标系转换方法

函数	说明
transform.wgs2gcj(lon, lat)	wgs84 → gcj-02
transform.gcj2wgs(lon, lat)	gcj-02 → wgs84
transform.gcj2bd(lon, lat)	gcj-02 → bd-09
transform.bd2gcj(lon, lat)	bd-09 → gcj-02
transform.bd2wgs(lon, lat)	bd-09 → wgs84
transform.wgs2bd(lon, lat)	wgs84 → bd-09

三、方法对比总结

方法	精度	适用场景	依赖库
geopy.geodesic	毫米级	全球范围，高精度	geopy
geopy.great_circle	中等	快速球面估算	geopy
haversine手动实现	中等	无库依赖	无
pyproj.geod	毫米级	复杂椭球模型计算	pyproj
平面近似	低	极小范围（城市内）	无