根据纬度/经度获取两点之间的距离-IT科技

摘要：问题描述：我尝试在基于纬度和经度查找距离中实现公式。该小程序对我测试的两点非常有用：但我的代码不起作用。from math import sin, cos, sqrt, atan2 R = 6373.0 lat1 = 52.2296756 lon1 = 21.0122287 lat2 = 52.4063...

问题描述：

我尝试在基于纬度和经度查找距离中实现公式。该小程序对我测试的两点非常有用：

在此处输入图片描述

但我的代码不起作用。

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。为什么？

解决方案 1：

自 GeoPy 版本 1.13 起， Vincenty 距离已被弃用-您应该改用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)

279.352901604将使用默认的椭球 WGS-84打印公里距离。（您也可以选择.miles或其他几种距离单位之一。）

解决方案 2：

值得注意的是，如果您只是需要一种快速简便的方法来找到两点之间的距离，我强烈建议您使用下面Kurt 的回答中描述的方法，而不是重新实现 Haversine——请参阅他的帖子了解基本原理。

这个答案仅侧重于回答 OP 遇到的特定错误。

这是因为在 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公里。

解决方案 3：

对于通过搜索引擎来到这里，并且只是在寻找一个开箱即用的解决方案的人（比如我），我建议安装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) # (latitude, longitude)
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) # (latitude, longitude)
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])

解决方案 4：

我找到了一个更简单、更强大的解决方案，即使用geodesic包geopy，因为无论如何您都很可能会在您的项目中使用它，所以不需要安装额外的包。

以下是我的解决方案：

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

地质学

解决方案 5：

有多种方法可以根据坐标（即纬度和经度）计算距离

安装并导入

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(radians(lat1)) * cos(radians(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 中使用 haversine

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

解决方案 6：

您可以使用Uber 的 H3函数point_dist()来计算两个（纬度、经度）点之间的球面距离。我们可以设置返回单位（“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

解决方案 7：

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)

解决方案 8：

在 2022 年，人们可以使用更新的 Python 库（即）发布解决此问题的混合 JavaScript 和 Python 代码geographiclib。总体好处是用户可以在现代设备上运行的网页上运行并查看结果。

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>

运行代码片段Hide results展开片段

解决方案 9：

（2022 年，实时 JavaScript 版本。）以下是使用较新的 JavaScript 库解决此问题的代码。总体好处是用户可以在现代设备上运行的网页上运行并查看结果。

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

运行代码片段Hide results展开片段

解决方案 10：

由于我们处于地理区域（！），因此这里有一个使用cartopy.geodesic.Geodesic.inverse()的解决方案。我相信您也可以使用pyproj.Geod.inv()来实现类似的效果。

注意：与此处提到的大多数其他包相反，cartopy 使用经纬度对。

import cartopy.geodesic as cgeo

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

globe = cgeo.Geodesic()

# Distance between two points
inv = globe.inverse(coords_1[::-1], coords_2[::-1])
print(inv.T[0]/1000)
# [279.3529016]

# Distances from one point to a list of other points
inv_2 = globe.inverse(
    coords_1[::-1],
    [coords_2[::-1], coords_3[::-1]],
)
print(inv_2.T[0]/1000)
# [279.3529016  750.45799898]

解决方案 11：

最简单的方法是使用haversine包。

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')

解决方案 12：

另一个有趣的用法是通过Pyodide和WebAssembly实现混合使用 JavaScript 和 Python，使用 Python 的库Pandas和geographiclib来获得解决方案也是可行的。

我花了额外的精力使用 Pandas 来准备输入数据，并在输出可用时将它们附加到solution列中。Pandas 为常见需求提供了许多有用的输入/输出功能。它的方法toHtml很方便在网页上呈现最终解决方案。

我发现此答案中的代码在某些iPhone和iPad设备上无法成功执行。但在较新的中端 Android 设备上它可以正常运行。

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

运行代码片段Hide results展开片段

解决方案 13：

这是我使用的东西，感谢Sviatoslav Oleksiv。

   earth_radius = {"km": 6371.0087714, "mile": 3959}

   earth_radius['km'] *
   acos(cos((math.radians(p1['latitude']))) *
   cos(math.radians(p2['latitude'])) *
   cos(math.radians(p2['longitude']) -
   math.radians(p1['longitude'])) +
   sin(math.radians(p1['latitude'])) *
   sin(math.radians(p2['latitude'])))