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

目录

• 一、基本方法
• 1. 使用geopy库（推荐）
• 2. 手动实现 Haversine 公式
• 3. 使用pyproj进行投影坐标系计算
• 4. 平面近似法（小范围适用）
• 二、不同坐标系的影响
• 常见坐标系
• coord-convert坐标系转换
• 三、方法对比总结
• 四、如何选择？
• 五、总结

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

以下是对应方法的示例。

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% 左右）。

适合无需外部库的场景。

示例代码

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结果一致。

适合需要平面坐标的高精度场景（如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一致。

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

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基础上二次加密的坐标系。

安装

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或无库平面近似。

以上为个人经验，希望能给大家一个参考，也希望大家多多支持风君子博客。

您可能感兴趣的文章:

• Python实现批量获取地址的经纬度
• python如何可视化经纬度信息
• Python实现读取.nc数据并提取指定时间与经纬度维度对应的变量数值
• Python的经纬度与xy坐标系相互转换方式
• python逐像素获取栅格经纬度分别保存在两个矩阵中(代码收藏)