我尝试在基于经纬度查找距离中实现公式。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。为什么?
我找到了一个更简单、更健壮的解决方案,即使用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
地珞
我找到了一个更简单、更健壮的解决方案,即使用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
地珞
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>
(2022年,JavaScript版本)下面是使用最新的JavaScript库解决这个问题的代码。总的好处是,用户可以在运行在现代设备上的web页面上看到结果。
// Using the 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` computation
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>