根据经纬度获取两点之间的距离

我尝试实现这个公式:http://andrew.hedges.name/experiments/haversine/ 这个小片对我正在测试的两点有好处:

enter image description here

但是我的代码没有工作。

from math import sin, cos, sqrt, atan2


R = 6373.0


lat1 = 52.2296756
lon1 = 21.0122287
lat2 = 52.406374
lon2 = 16.9251681


dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2
c = 2 * atan2(sqrt(a), sqrt(1-a))
distance = R * c


print "Result", distance
print "Should be", 278.546

它返回的距离是5447.05546147。为什么?

435030 次浏览
小开
最佳答案

编辑:只是一个注释，如果你只是需要一个快速简单的方法来找到两点之间的距离，我强烈建议使用下面库尔特的回答中描述的方法，而不是重新实现Haversine——查看他的帖子来了解基本原理。

这个答案只关注于回答OP遇到的特定错误。

这是因为在Python中，所有的三角函数使用弧度，而不是度数。

你可以手动将数字转换为弧度，或者使用math模块中的radians函数:

from math import sin, cos, sqrt, atan2, radians


# approximate radius of earth in km
R = 6373.0


lat1 = radians(52.2296756)
lon1 = radians(21.0122287)
lat2 = radians(52.406374)
lon2 = radians(16.9251681)


dlon = lon2 - lon1
dlat = lat2 - lat1


a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))


distance = R * c


print("Result:", distance)
print("Should be:", 278.546, "km")

距离现在返回正确的值278.545589351 km。

小开

对于像我这样通过搜索引擎来这里并只是寻找一个开箱即用的解决方案的人，我建议安装mpu。通过pip install mpu --user安装它，并像这样使用它来获得半正矢距离:

import mpu


# Point one
lat1 = 52.2296756
lon1 = 21.0122287


# Point two
lat2 = 52.406374
lon2 = 16.9251681


# What you were looking for
dist = mpu.haversine_distance((lat1, lon1), (lat2, lon2))
print(dist)  # gives 278.45817507541943.

另一个包是gpxpy。

如果你不想要依赖，你可以使用:

import math




def distance(origin, destination):
"""
Calculate the Haversine distance.


Parameters
----------
origin : tuple of float
(lat, long)
destination : tuple of float
(lat, long)


Returns
-------
distance_in_km : float


Examples
--------
>>> origin = (48.1372, 11.5756)  # Munich
>>> destination = (52.5186, 13.4083)  # Berlin
>>> round(distance(origin, destination), 1)
504.2
"""
lat1, lon1 = origin
lat2, lon2 = destination
radius = 6371  # km


dlat = math.radians(lat2 - lat1)
dlon = math.radians(lon2 - lon1)
a = (math.sin(dlat / 2) * math.sin(dlat / 2) +
math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *
math.sin(dlon / 2) * math.sin(dlon / 2))
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
d = radius * c


return d




if __name__ == '__main__':
import doctest
doctest.testmod()

另一个可选包是haversine

from haversine import haversine, Unit


lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)


haversine(lyon, paris)
>> 392.2172595594006  # in kilometers


haversine(lyon, paris, unit=Unit.MILES)
>> 243.71201856934454  # in miles


# you can also use the string abbreviation for units:
haversine(lyon, paris, unit='mi')
>> 243.71201856934454  # in miles


haversine(lyon, paris, unit=Unit.NAUTICAL_MILES)
>> 211.78037755311516  # in nautical miles

他们声称对两个向量中所有点之间的距离进行了性能优化

from haversine import haversine_vector, Unit


lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)
new_york = (40.7033962, -74.2351462)


haversine_vector([lyon, lyon], [paris, new_york], Unit.KILOMETERS)


>> array([ 392.21725956, 6163.43638211])
小开

更新:04/2018: Vincenty distance改为自GeoPy 1.13版后已弃用 - 你应该使用geopy.distance.distance() !

上面的答案是基于半正矢公式，它假设地球是一个球体，这导致误差高达0.5%(根据help(geopy.distance))。Vincenty距离使用更精确的椭球模型，如wgs - 84，并在geopy中实现。例如,

import geopy.distance


coords_1 = (52.2296756, 21.0122287)
coords_2 = (52.406374, 16.9251681)


print geopy.distance.geodesic(coords_1, coords_2).km

将使用默认椭球WGS-84打印279.352901604千米的距离。(你也可以选择.miles或其他距离单位之一)。

小开 
import numpy as np




def Haversine(lat1,lon1,lat2,lon2, **kwarg):
"""
This uses the ‘haversine’ formula to calculate the great-circle distance between two points – that is,
the shortest distance over the earth’s surface – giving an ‘as-the-crow-flies’ distance between the points
(ignoring any hills they fly over, of course!).
Haversine
formula:    a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
c = 2 ⋅ atan2( √a, √(1−a) )
d = R ⋅ c
where   φ is latitude, λ is longitude, R is earth’s radius (mean radius = 6,371km);
note that angles need to be in radians to pass to trig functions!
"""
R = 6371.0088
lat1,lon1,lat2,lon2 = map(np.radians, [lat1,lon1,lat2,lon2])


dlat = lat2 - lat1
dlon = lon2 - lon1
a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2) **2
c = 2 * np.arctan2(a**0.5, (1-a)**0.5)
d = R * c
return round(d,4)
小开

我得到了一个更简单和健壮的解决方案，这是使用geopy包中的geodesic，因为你很可能在你的项目中使用它，所以不需要额外的包安装。

以下是我的解决方案:

from geopy.distance import geodesic




origin = (30.172705, 31.526725)  # (latitude, longitude) don't confuse
dist = (30.288281, 31.732326)


print(geodesic(origin, dist).meters)  # 23576.805481751613
print(geodesic(origin, dist).kilometers)  # 23.576805481751613
print(geodesic(origin, dist).miles)  # 14.64994773134371

geopy

小开

有多种方法来计算基于坐标的距离，即纬度和经度

安装和导入

from geopy import distance
from math import sin, cos, sqrt, atan2, radians
from sklearn.neighbors import DistanceMetric
import osrm
import numpy as np

定义坐标

lat1, lon1, lat2, lon2, R = 20.9467,72.9520, 21.1702, 72.8311, 6373.0
coordinates_from = [lat1, lon1]
coordinates_to = [lat2, lon2]

使用半正矢

dlon = radians(lon2) - radians(lon1)
dlat = radians(lat2) - radians(lat1)
    

a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
    

distance_haversine_formula = R * c
print('distance using haversine formula: ', distance_haversine_formula)

使用哈弗辛和sklearn

dist = DistanceMetric.get_metric('haversine')
    

X = [[radians(lat1), radians(lon1)], [radians(lat2), radians(lon2)]]
distance_sklearn = R * dist.pairwise(X)
print('distance using sklearn: ', np.array(distance_sklearn).item(1))

使用OSRM

osrm_client = osrm.Client(host='http://router.project-osrm.org')
coordinates_osrm = [[lon1, lat1], [lon2, lat2]] # note that order is lon, lat
    

osrm_response = osrm_client.route(coordinates=coordinates_osrm, overview=osrm.overview.full)
dist_osrm = osrm_response.get('routes')[0].get('distance')/1000 # in km
print('distance using OSRM: ', dist_osrm)

使用geopy

distance_geopy = distance.distance(coordinates_from, coordinates_to).km
print('distance using geopy: ', distance_geopy)
    

distance_geopy_great_circle = distance.great_circle(coordinates_from, coordinates_to).km
print('distance using geopy great circle: ', distance_geopy_great_circle)

输出

distance using haversine formula:  26.07547017310917
distance using sklearn:  27.847882224769783
distance using OSRM:  33.091699999999996
distance using geopy:  27.7528030550408
distance using geopy great circle:  27.839182219511834
小开

你可以使用超级的H3，point_dist()函数来计算两个(lat, lng)点之间的球面距离。我们可以设置返回单位('km'， 'm'，或'rads')。默认单位为Km。

例子:

import h3


coords_1 = (52.2296756, 21.0122287)
coords_2 = (52.406374, 16.9251681)
distance = h3.point_dist(coords_1, coords_2, unit='m') # to get distance in meters

希望这对你有用!

小开

最简单的方法是用哈弗辛包装。


import haversine as hs




coord_1 = (lat, lon)
coord_2 = (lat, lon)
x = hs.haversine(coord_1,coord_2)
print(f'The distance is {x} km')
小开

(2022年，实时javascript版本)下面是使用最新的javascript库解决这个问题的代码。总的好处是，用户可以在运行在现代设备上的web页面上看到结果。

// Using WGS84 ellipsoid model for computation
var geod84 = geodesic.Geodesic.WGS84;
// Input data
lat1 = 52.2296756;
lon1 = 21.0122287;
lat2 = 52.406374;
lon2 = 16.9251681;
// Do the classic `geodetic inversion` computatioin
geod84inv = geod84.Inverse(lat1, lon1, lat2, lon2);
// Present the solution (only the geodetic distance)
console.log("The distance is " + (geod84inv.s12/1000).toFixed(5) + " km.");
<script type="text/javascript" src="https://cdn.jsdelivr.net/npm/geographiclib-geodesic@2.0.0/geographiclib-geodesic.min.js">
</script>

小开

在2022年，人们可以发布混合javascript+python代码，使用最新的python库，即geographiclib来解决这个问题。总的好处是，用户可以在运行在现代设备上的web页面上看到结果。

async function main(){
let pyodide = await loadPyodide();
await pyodide.loadPackage(["micropip"]);
  

console.log(pyodide.runPythonAsync(`
import micropip
await micropip.install('geographiclib')
from geographiclib.geodesic import Geodesic
lat1 = 52.2296756;
lon1 = 21.0122287;
lat2 = 52.406374;
lon2 = 16.9251681;
ans = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2)
dkm = ans["s12"] / 1000
print("Geodesic solution", ans)
print(f"Distance = {dkm:.4f} km.")
`));
  

}


main();
<script src="https://cdn.jsdelivr.net/pyodide/v0.21.0/full/pyodide.js"></script>

小开
通过pyodide和webassembly实现混合javascript+python来使用Python的库pandas+geographiclib获得解决方案的另一种有趣的用法也是可行的。 我额外地使用pandas来准备输入数据，当输出可用时，将它们附加到solution列。Pandas为常见需求提供了许多有用的输入/输出特性。它的方法toHtml可以方便地在web页面

上显示最终解决方案

< em >编辑 我发现这个答案中的代码在某些iphone和ipad设备上执行不成功。但在更新的中端Android设备上运行也很好。我会找到一个方法来纠正这个问题，并尽快更新

我的旁注，我知道我的答案不像其他答案那样直接回答OP问题。但最近外界表示，stackoverflow中的许多答案已经过时，并试图引导人们远离这里。

async function main(){
let pyodide = await loadPyodide();
await pyodide.loadPackage(["pandas", "micropip"]);
console.log(pyodide.runPythonAsync(`
import micropip
import pandas as pd
import js
print("Pandas version: " + pd.__version__)
await micropip.install('geographiclib')
from geographiclib.geodesic import Geodesic
import geographiclib as gl
print("Geographiclib version: " + gl.__version__)
data = {'Description': ['Answer to the question', 'Bangkok to Tokyo'],
'From_long': [21.0122287, 100.6],
'From_lat': [52.2296756, 13.8],
'To_long': [16.9251681, 139.76],
'To_lat': [52.406374, 35.69],
'Distance_km': [0, 0]}
df1 = pd.DataFrame(data)
collist = ['Description','From_long','From_lat','To_long','To_lat']
div2 = js.document.createElement("div")
div2content = df1.to_html(buf=None, columns=collist, col_space=None, header=True, index=True)
div2.innerHTML = div2content
js.document.body.append(div2)
arr="<i>by Swatchai</i>"


def dkm(frLat,frLon,toLat,toLon):
print("frLon,frLat,toLon,toLat:", frLon, "|", frLat, "|", toLon, "|", toLat)
dist = Geodesic.WGS84.Inverse(frLat, frLon, toLat, toLon)
return dist["s12"] / 1000
  

collist = ['Description','From_long','From_lat','To_long','To_lat','Distance_km']
dist = []
for ea in zip(df1['From_lat'].values, df1['From_long'].values, df1['To_lat'].values, df1['To_long'].values):
ans = dkm(*ea)
print("ans=", ans)
dist.append(ans)
  

df1['Distance_km'] = dist
# Update content
div2content = df1.to_html(buf=None, columns=collist, col_space=None, header=True, index=False)
div2.innerHTML = div2content
js.document.body.append(div2)


# Using Haversine Formula
from math import sin, cos, sqrt, atan2, radians, asin
# approximate radius of earth in km from wikipedia
R = 6371
lat1 = radians(52.2296756)
lon1 = radians(21.0122287)
lat2 = radians(52.406374)
lon2 = radians(16.9251681)
# https://en.wikipedia.org/wiki/Haversine_formula
def hav(angrad):
return (1-cos(angrad))/2
h = hav(lat2-lat1)+cos(lat2)*cos(lat1)*hav(lon2-lon1)
dist2 = 2*R*asin(sqrt(h))
print(f"Distance by haversine formula = {dist2:8.6f} km.")
`));




}
main();
<script src="https://cdn.jsdelivr.net/pyodide/v0.21.0/full/pyodide.js"></script>
Pyodide implementation<br>


提交答案