Python中经纬度距离计算的实现方式
作者：AI手记叨叨
文章介绍Python中计算经纬度距离的方法及中国加密坐标系转换工具,主要方法包括geopy（Vincenty/Karney）、Haversine、pyproj和平面近似法,适用于不同精度和场景需求,需注意坐标系转换（如GCJ-02→WGS84）的不可逆性及误差问题
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 |
| 平面近似 | 低 | 极小范围（城市内） | 无 |
四、如何选择？
- 通用场景：直接使用
geopy.geodesic（平衡精度和易用性）。 - 学术与工程：需要椭球模型时用
pyproj。 - 轻量级需求：手动Haversine或无库平面近似。
五、总结
以上为个人经验，希望能给大家一个参考，也希望大家多多支持脚本之家。