我尝试在基于经纬度查找距离中实现公式。applet对我正在测试的两点很好:
但是我的代码没有工作。
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。为什么?
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)
在2022年,人们可以发布JavaScript和Python混合代码,使用最新的Python库,即地理库来解决这个问题。总的好处是,用户可以在运行在现代设备上的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>
我找到了一个更简单、更健壮的解决方案,即使用geogeoy包中的测地线,因为你很可能在你的项目中使用它,所以不需要额外的包安装。
以下是我的解决方案:
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
地珞
注意,如果你只是需要一种快速简单的方法来找到两点之间的距离,我强烈建议使用Kurt回答中描述的方法,而不是重新实现haversine,查看他的帖子来了解基本原理。
这个答案只关注OP遇到的特定bug。
这是因为在Python中,所有的三角函数都使用弧度,而不是度。
您可以手动将数字转换为弧度,或使用数学模块中的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。