根据经纬度坐标计算两点距离

博主：稻叶
发布时间：2023 年 10 月 22 日
232 次浏览
暂无评论
3242字数
分类：应用

# python 实现 haversine 公式计算

根据经纬度计算两点距离，可以利用半正矢公式来进行计算。

## 半正矢公式

半正矢公式是一种根据两点的经度和纬度来确定 大圆距离[^1] 的计算方法，是导航学中的一种基本方法。该公式使用正半矢函数表达，由 $haversin(\theta)=sin^2(\frac{\theta}{2})$ 而来。

设球面上任意两点之间的圆心角 $\theta$:

$\theta=\frac{d}{r}$

* $d$ :两点之间的距离
* $r$ :地球的半径

$$
hav(\theta)=hav(\varphi_2-\varphi_1)+cos(\varphi_1)cos(\varphi_2)hav(\lambda_2-\lambda_1)
$$

公式中：

* $\varphi$ :纬度(*lat.*)
* $\lambda$ :经度(*long.*)
* $hav()$ :半正矢函数(*haversin*)

$$
  hav(\theta)=sin^2(\frac{\theta}{2})=\frac{1-cos(\theta)}{2}
  $$

变形得到：

$$
d=2r\,arcsin\sqrt{sin^2(\frac{\varphi_2-\varphi_1}{2})+cos(\varphi_1)cos(\varphi_2)sin^2(\frac{\lambda_2-\lambda_1}{2})}
$$

之后我们便可以利用该公式进行地理空间距离计算了

## 实现

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

EARTH_RADIUS = 6371  # 地球平均半径：6371km

def get_distance_hav(lon0, lat0, lon1, lat1):
    # 用haversine公式计算球面两点间的距离
    lat0 = radians(lat0)
    lat1 = radians(lat1)
    lon0 = radians(lon0)
    lon1 = radians(lon1)

dlat = fabs(lat1 - lat0)
    dlon = fabs(lon1 - lon0)
    h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlon)
    distance = 2 * EARTH_RADIUS * asin(sqrt(h))

return distance # 单位：km

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

# 利用 GeoPy 库计算

Geopy 是一个 Python 库，主要用于处理地理空间数据和地理位置信息。它提供了一系列功能，用于地理编码（地址到经纬度的转换）、反向地理编码（经纬度到地址的转换）、计算两点之间的距离等地理信息处理任务。

geopy 的距离计算功能在`geopy.distance`中。其计算距离有两种方法：大圆距离([great-circle distance](https://en.wikipedia.org/wiki/Great-circle_distance))与测地线距离([geodesic distance](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid))。

其区别在于：

* 大圆距离将地球简化为一个半径为`6371.009km`的圆球，其计算的距离是球面上两点间的大圆距离；
* 测地线距离使用目前国际通用的方法，用旋转椭球面表示地球，其计算的距离是两点在椭球面上的最短距离。

`geopy.distance.distance`默认使用更精确的`geodesic`。

## 安装

```shell
pip install geopy
```

## 计算

`distance`以元组`(lat, lon)`传入经纬度

```python
from geopy import distance

# 定义两个坐标点的经纬度
coord0 = (lat0, lon0)
coord1 = (lat1, lon1)

# distance
result = distance.distance(coord0,coord1).km
# geodesic
result1 = distance.geodesic(coord0,coord1).km
# 更改使用的椭球体模型为"Airy (1830)"
result2 = distance.geodesic(coord0,coord1,ellipsoid="Airy (1830)").km
# 使用 great_circle 距离
result3 = distance.great_circle(coord0,coord1).km
```

geodesic有许多流行的椭球体模型，哪一个模型最精确取决于你的点在地球上的位置。默认值是WGS-84椭球体，它是全局最精确的。

使用模型名称从`ELLIPSOIDS`字典选取值：

```python
ELLIPSOIDS = {
    # model           major (km)   minor (km)     flattening
    'WGS-84':        (6378.137, 6356.7523142, 1 / 298.257223563),
    'GRS-80':        (6378.137, 6356.7523141, 1 / 298.257222101),
    'Airy (1830)':   (6377.563396, 6356.256909, 1 / 299.3249646),
    'Intl 1924':     (6378.388, 6356.911946, 1 / 297.0),
    'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465),
    'GRS-67':        (6378.1600, 6356.774719, 1 / 298.25)
}
```

更多关于GeoPy的使用可以查看[GeoPy文档](https://www.osgeo.cn/geopy/)

# 计算结果

选取两个坐标点进行计算，使用`%timeit`计算函数运行时间

```python
lon0, lat0, lon1, lat1 = (113.324548, 23.106413, 113.322552, 23.115056)
coord0 = (lat0, lon0)
coord1 = (lat1, lon1)

result = get_distance_hav(lon0,lat0,lon1,lat1)
result1 = distance.geodesic(coord0,coord1).km
result2 = distance.great_circle(coord0,coord1).km

print("{:<15}{:.6f} km".format("haversine:",result))
print("{:<15}{:.6f} km".format("geodesic:",result1))
print("{:<15}{:.6f} km".format("great_circle:",result2))
```

输出结果：

```shell
haversine:     0.9825 km
geodesic:      0.9788 km
great_circle:  0.9825 km
```

平均耗时：

```shell
haversine:     1.31 µs
geodesic:      131 µs
great_circle:  8.62 µs
```

根据输出结果可以看出，haversine与geodesic的计算结果的差异较小，但前者的计算时间更短，在不需要高度精确的场景下，牺牲一定精度可以换取更高的效率。

[^1]: 指从球面的一点A出发到达球面上另一点B，所经过的最短路径的长度。

最后修改：2023 年 10 月 23 日

如果觉得我的文章对你有用，请随意赞赏

发表评论取消回复
使用cookie技术保留您的个人信息以便您下次快速评论，继续评论表示您已同意该条款

评论 *

私密评论

名称 *

🎲

邮箱 *

地址

根据经纬度坐标计算两点距离

稻叶 • 2023 年 10 月 22 日

# python 实现 haversine 公式计算

根据经纬度计算两点距离，可以利用半正矢公式来进行计算。

## 半正矢公式

设球面上任意两点之间的圆心角 $\theta$:

$\theta=\frac{d}{r}$

* $d$ :两点之间的距离
* $r$ :地球的半径

$$
hav(\theta)=hav(\varphi_2-\varphi_1)+cos(\varphi_1)cos(\varphi_2)hav(\lambda_2-\lambda_1)
$$

公式中：

* $\varphi$ :纬度(*lat.*)
* $\lambda$ :经度(*long.*)
* $hav()$ :半正矢函数(*haversin*)

$$
  hav(\theta)=sin^2(\frac{\theta}{2})=\frac{1-cos(\theta)}{2}
  $$

变形得到：

$$
d=2r\,arcsin\sqrt{sin^2(\frac{\varphi_2-\varphi_1}{2})+cos(\varphi_1)cos(\varphi_2)sin^2(\frac{\lambda_2-\lambda_1}{2})}
$$

之后我们便可以利用该公式进行地理空间距离计算了

## 实现

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

EARTH_RADIUS = 6371  # 地球平均半径：6371km

dlat = fabs(lat1 - lat0)
    dlon = fabs(lon1 - lon0)
    h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlon)
    distance = 2 * EARTH_RADIUS * asin(sqrt(h))

return distance # 单位：km

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

# 利用 GeoPy 库计算

其区别在于：

`geopy.distance.distance`默认使用更精确的`geodesic`。

## 安装

```shell
pip install geopy
```

## 计算

`distance`以元组`(lat, lon)`传入经纬度

```python
from geopy import distance

# 定义两个坐标点的经纬度
coord0 = (lat0, lon0)
coord1 = (lat1, lon1)

geodesic有许多流行的椭球体模型，哪一个模型最精确取决于你的点在地球上的位置。默认值是WGS-84椭球体，它是全局最精确的。

使用模型名称从`ELLIPSOIDS`字典选取值：

更多关于GeoPy的使用可以查看[GeoPy文档](https://www.osgeo.cn/geopy/)

# 计算结果

选取两个坐标点进行计算，使用`%timeit`计算函数运行时间

```python
lon0, lat0, lon1, lat1 = (113.324548, 23.106413, 113.322552, 23.115056)
coord0 = (lat0, lon0)
coord1 = (lat1, lon1)

result = get_distance_hav(lon0,lat0,lon1,lat1)
result1 = distance.geodesic(coord0,coord1).km
result2 = distance.great_circle(coord0,coord1).km

print("{:<15}{:.6f} km".format("haversine:",result))
print("{:<15}{:.6f} km".format("geodesic:",result1))
print("{:<15}{:.6f} km".format("great_circle:",result2))
```

输出结果：

```shell
haversine:     0.9825 km
geodesic:      0.9788 km
great_circle:  0.9825 km
```

平均耗时：

```shell
haversine:     1.31 µs
geodesic:      131 µs
great_circle:  8.62 µs
```

[^1]: 指从球面的一点A出发到达球面上另一点B，所经过的最短路径的长度。

发表评论 取消回复 使用cookie技术保留您的个人信息以便您下次快速评论，继续评论表示您已同意该条款

根据经纬度坐标计算两点距离

发表评论取消回复
使用cookie技术保留您的个人信息以便您下次快速评论，继续评论表示您已同意该条款