Python 优雅地利用两点经纬度计算地理空间距离

主要内容概览：

 

一、根据经纬度计算距离的基本原理

处理地理数据时，经常需要用到两个地理位置间的距离。以点（经度值，纬度值）表示地点的经纬度坐标，比如 A 点经纬度（30.553949，114.357399），B点经纬度（129.1344，25.5465），求 AB 两点之间的距离。

已知地球上任意两点（lng1，lat1），（lng2， lat2）的经纬度坐标，求两点间的距离可以利用 haversine 公式：

• 首先，将经纬度坐标的角度化成弧度（rlng1，rlat1），（rlng2，rlat2)
• 然后，利用如下公式计算：

其中， a 表示两点纬度的差值，即 a = rlat1 - rlat2； b表示两点经度的差值，即 b = rlng1 - rlng2，其中 r 表示地球的半径。

以计算清华大学与北京大学距离为例，利用百度地图（ 网址：百度地图）测距可知，清华大学与北京大学距离大约为1.8km，如下所示：

二、获取地点的经纬度

利用高德地图地理编码 / 逆地理编码 API ，来获取北京大学和清华大学的经纬度，Python代码实现如下：

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import requests


KEY = '你的高德地图应用访问码'


def geocode(address):
    """
    根据地名，获取经纬度
    """
    parameters = {'address': address, 'key': KEY}
    base = 'http://restapi.amap.com/v3/geocode/geo'
    response = requests.get(base, parameters)
    answer = response.json()
    # print(answer)
    # print(type(answer))
    print(address + "的经纬度：", answer['geocodes'][0]['location'])


if __name__ == "__main__":
    address1 = '北京大学'
    lon_lat1 = geocode(address1)
    address2 = '清华大学'
    lon_lat2 = geocode(address2)

 其中，KEY为你自己申请的高德地图应用授权访问码。只需输入地址参数，就可以打印出对应地址的经纬度坐标值。

运行结果为：

北京大学的经纬度： 116.308264,39.995304
清华大学的经纬度： 116.326759,40.003304

三、手写 haversine 公式计算

根据 haversine 公式，自定义实现python版本计算公式，如下：

from math import sin, asin, cos, radians, fabs, sqrt

EARTH_RADIUS = 6371.137      # 地球平均半径大约6371km


def hav(theta):
    s = sin(theta / 2)
    return s * s


def get_distance_hav(lat0, lng0, lat1, lng1):
    # 用haversine公式计算球面两点间的距离
    # 经纬度转换成弧度
    lat0 = radians(lat0)
    lat1 = radians(lat1)
    lng0 = radians(lng0)
    lng1 = radians(lng1)
    dlng = fabs(lng0 - lng1)
    dlat = fabs(lat0 - lat1)
    h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlng)
    distance = 2 * EARTH_RADIUS * asin(sqrt(h))      # km
    return distance

if __name__ == "__main__":
    result = get_distance_hav(39.995304, 116.308264, 40.003304, 116.326759)
    print("距离：{:.3}km".format(result))

运行结果为：

距离：1.809km

四、利用 geopy 库计算

geopy是一个关于地理编码的 Python 库。主要有以下几个功能：

• 地理编码：将字符串转换为地理位置。
• 逆地理编码：用于将地理坐标转换为具体地址。
• 计算两个点的距离：经纬度距离和球面距离。

pip install安装geopy库，命令为：

pip install geopy -i http://pypi.douban.com/simple --trusted-host pypi.douban.com

Python实现代码如下所示：

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from geopy.distance import geodesic


if __name__ == "__main__":
    distance = geodesic((39.995304, 116.308264), (40.003304, 116.326759)).km  
    print("距离：{:.3f}km".format(distance))

运行结果如下所示：

距离：1.812km

手写 haversine 公式计算距离为：1.809km，而直接调用 geopy 库的 geodesic 方法计算距离为：1.812km，二者计算结果相差很小。 

更多有关 geopy 的使用，可以学习官方文档：欢迎使用Geopy文档！ — GeoPy 2.2.0 文档。

五、利用 haversine库 计算

pypi地址：haversine · PyPI

Calculate the distance (in various units) between two points on Earth using their latitude and longitude.

用经纬度计算地球上两点之间的距离（以不同单位表示），pip install安装即可，命令为：

pip install haversine -i http://pypi.douban.com/simple --trusted-host pypi.douban.com

Python实现代码如下：

from haversine import haversine, Unit


if __name__ == "__main__":
    # 两点的经纬度
    point1 = (39.995304, 116.308264)
    point2 = (40.003304, 116.326759)
    result1 = haversine(point1, point2, unit=Unit.KILOMETERS)    # km
    result2 = haversine(point1, point2, unit=Unit.METERS)        # m
    # 打印计算结果
    print("距离：{:.3f}km".format(result1))
    print("距离：{:.3f}m".format(result2))

运行结果为：

距离：1.809km
距离：1809.223m

可以发现，利用 haversine 库计算距离的结果与手写 haversine 公式计算的结果一致。