我尝试在基于经纬度查找距离中实现公式。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。为什么?


当前回答

(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>

其他回答

在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>

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

安装和导入

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

注意,如果你只是需要一种快速简单的方法来找到两点之间的距离,我强烈建议使用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。

Vincenty距离从GeoPy 1.13版开始就被弃用了-你应该使用geo .distance.distance()来代替!


上面的答案是基于haversine公式,该公式假设地球是一个球体,结果误差高达0.5%(根据help(earth .distance))。Vincenty距离采用更精确的椭球模型,如WGS-84,并在地质学中实现。例如,

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或其他距离单位。)

(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>