计算两个经纬度点之间的距离？ （Haversine 公式）
Posted
技术标签:
【中文标题】计算两个经纬度点之间的距离？ （Haversine 公式）【英文标题】：Calculate distance between two latitude-longitude points? (Haversine formula) 【发布时间】：2010-09-06 21:00:41 【问题描述】：如何计算经纬度指定的两点之间的距离？
为了澄清，我想以公里为单位；这些点使用 WGS84 系统，我想了解可用方法的相对准确性。
【问题讨论】：
为了更准确 - 请参阅***.com/questions/1420045/… 请注意，您不能在像 WGS 84 这样的旋转椭球体上应用半正弦公式。您只能在具有半径的球体上应用此方法。 这里的大多数答案都是使用简单的球面三角法，因此与 GPS 系统中使用的 WGS84 椭球距离相比，结果相当粗糙。一些答案确实参考了 Vincenty 的椭球公式，但该算法是为 1960 年代的台式计算器设计的，并且存在稳定性和准确性问题；我们现在有更好的硬件和软件。请参阅GeographicLib 以了解具有多种语言实现的高质量库。 @MikeT - 虽然这里的许多答案似乎在小距离内很有用：如果您从 WGS 84 获取 lat/long，并应用 Haversine 就好像那些 是球体上的点，你没有得到答案，其错误只是由于地球的扁平化因素，所以可能在更准确公式的 1% 以内？需要注意的是，这些距离很短，比如在一个城镇内。 对于这些平台：Mono/.NET 4.5/.NET Core/Windows Phone 8.x/通用 Windows 平台/Xamarin ios/Xamarin android 参见***.com/a/54296314/2736742 【参考方案1】：这是另一个转换为 Ruby 的代码：
include Math
#Note: from/to = [lat, long]
def get_distance_in_km(from, to)
radians = lambda |deg| deg * Math.PI / 180
radius = 6371 # Radius of the earth in kilometer
dLat = radians[to[0]-from[0]]
dLon = radians[to[1]-from[1]]
cosines_product = Math.sin(dLat/2) * Math.sin(dLat/2) + Math.cos(radians[from[0]]) * Math.cos(radians[to[1]]) * Math.sin(dLon/2) * Math.sin(dLon/2)
c = 2 * Math.atan2(Math.sqrt(cosines_product), Math.sqrt(1-cosines_product))
return radius * c # Distance in kilometer
end
【讨论】：
【参考方案2】：由于这是该主题最受欢迎的讨论，我将在此处添加我从 2019 年末至 2020 年初的经验。添加到现有答案 - 我的重点是找到一个准确且快速（即矢量化）的解决方案。
让我们从这里的答案最常用的方法开始 - Haversine 方法。向量化很简单，见下面的python示例：
def haversine(lat1, lon1, lat2, lon2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
All args must be of equal length.
Distances are in meters.
Ref:
https://***.com/questions/29545704/fast-haversine-approximation-python-pandas
https://ipython.readthedocs.io/en/stable/interactive/magics.html
"""
Radius = 6.371e6
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
c = 2 * np.arcsin(np.sqrt(a))
s12 = Radius * c
# initial azimuth in degrees
y = np.sin(lon2-lon1) * np.cos(lat2)
x = np.cos(lat1)*np.sin(lat2) - np.sin(lat1)*np.cos(lat2)*np.cos(dlon)
azi1 = np.arctan2(y, x)*180./math.pi
return 's12':s12, 'azi1': azi1
在准确性方面，它是最不准确的。***指出平均 0.5% 的相对偏差没有任何来源。我的实验表明偏差较小。下面是 100,000 个随机点与我的库的比较，应该精确到毫米级：
np.random.seed(42)
lats1 = np.random.uniform(-90,90,100000)
lons1 = np.random.uniform(-180,180,100000)
lats2 = np.random.uniform(-90,90,100000)
lons2 = np.random.uniform(-180,180,100000)
r1 = inverse(lats1, lons1, lats2, lons2)
r2 = haversine(lats1, lons1, lats2, lons2)
print("Max absolute error: :4.2fm".format(np.max(r1['s12']-r2['s12'])))
print("Mean absolute error: :4.2fm".format(np.mean(r1['s12']-r2['s12'])))
print("Max relative error: :4.2f%".format(np.max((r2['s12']/r1['s12']-1)*100)))
print("Mean relative error: :4.2f%".format(np.mean((r2['s12']/r1['s12']-1)*100)))
输出：
Max absolute error: 26671.47m
Mean absolute error: -2499.84m
Max relative error: 0.55%
Mean relative error: -0.02%
因此，在 100,000 个随机坐标对上平均有 2.5 公里的偏差，这可能对大多数情况都有好处。
下一个选项是 Vincenty 的公式，该公式精确到毫米，具体取决于收敛标准，并且还可以矢量化。它确实存在在对跖点附近收敛的问题。您可以通过放宽收敛标准使其在这些点收敛，但准确度会下降到 0.25% 甚至更多。在对映点之外，Vincenty 将提供接近 Geographiclib 的结果，平均相对误差小于 1.e-6。
这里提到的 Geographiclib 确实是当前的黄金标准。它有几个实现并且相当快，尤其是如果您使用的是 C++ 版本。
现在，如果您打算将 Python 用于超过 10k 点的任何内容，我建议您考虑我的矢量化实现。我根据自己的需要创建了一个带有矢量化 Vincenty 例程的 geovectorslib 库，它使用 Geographiclib 作为近对跖点的后备。下面是 100k 点与 Geographiclib 的比较。如您所见，它为 100k 点提供了 20 倍的逆向改进和 100 倍的直接改进，并且差距会随着点数的增加而增长。精度方面，它将在 Georgraphiclib 的 1.e-5 rtol 内。
Direct method for 100,000 points
94.9 ms ± 25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
9.79 s ± 1.4 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
Inverse method for 100,000 points
1.5 s ± 504 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
24.2 s ± 3.91 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
【讨论】：
【参考方案3】：计算距离（尤其是大距离）的主要挑战之一是考虑地球的曲率。如果地球是平的，计算两点之间的距离就像计算直线一样简单！ Haversine 公式包括一个表示地球半径的常数（它是下面的 R 变量）。根据您是以英里还是公里为单位，它分别等于 3956 英里或 6367 公里。 基本公式为：
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 (where R is the radius of the Earth) R = 6367 km OR 3956 mi
lat1, lon1: The Latitude and Longitude of point 1 (in decimal degrees)
lat2, lon2: The Latitude and Longitude of point 2 (in decimal degrees)
unit: The unit of measurement in which to calculate the results where:
'M' is statute miles (default)
'K' is kilometers
'N' is nautical miles
样本
function distance(lat1, lon1, lat2, lon2, unit)
try
var radlat1 = Math.PI * lat1 / 180
var radlat2 = Math.PI * lat2 / 180
var theta = lon1 - lon2
var radtheta = Math.PI * theta / 180
var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
dist = Math.acos(dist)
dist = dist * 180 / Math.PI
dist = dist * 60 * 1.1515
if (unit == "K")
dist = dist * 1.609344
if (unit == "N")
dist = dist * 0.8684
return dist
catch (err)
console.log(err);
【讨论】：
虽然此代码可能会解决问题，including an explanation 关于如何以及为什么解决问题将真正有助于提高您的帖子质量，并可能导致更多的赞成票。请记住，您正在为将来的读者回答问题，而不仅仅是现在提问的人。请edit您的回答添加解释并说明适用的限制和假设。【参考方案4】：如果您想要行驶距离/路线（在此处发布，因为这是谷歌上两点之间距离的第一个结果，但对于大多数人来说行驶距离更有用），您可以使用Google Maps Distance Matrix Service：
getDrivingDistanceBetweenTwoLatLong(origin, destination)
return new Observable(subscriber =>
let service = new google.maps.DistanceMatrixService();
service.getDistanceMatrix(
origins: [new google.maps.LatLng(origin.lat, origin.long)],
destinations: [new google.maps.LatLng(destination.lat, destination.long)],
travelMode: 'DRIVING'
, (response, status) =>
if (status !== google.maps.DistanceMatrixStatus.OK)
console.log('Error:', status);
subscriber.error(error: status, status: status);
else
console.log(response);
try
let valueInMeters = response.rows[0].elements[0].distance.value;
let valueInKms = valueInMeters / 1000;
subscriber.next(valueInKms);
subscriber.complete();
catch(error)
subscriber.error(error: error, status: status);
);
);
【讨论】：
【参考方案5】：function getDistanceFromLatLonInKm(position1, position2)
"use strict";
var deg2rad = function (deg) return deg * (Math.PI / 180); ,
R = 6371,
dLat = deg2rad(position2.lat - position1.lat),
dLng = deg2rad(position2.lng - position1.lng),
a = Math.sin(dLat / 2) * Math.sin(dLat / 2)
+ Math.cos(deg2rad(position1.lat))
* Math.cos(deg2rad(position2.lat))
* Math.sin(dLng / 2) * Math.sin(dLng / 2),
c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return R * c;
console.log(getDistanceFromLatLonInKm(
lat: 48.7931459, lng: 1.9483572,
lat: 48.827167, lng: 2.2459745
));
【讨论】：
你打错了，第二个cos应该是pos2【参考方案6】：LUA 中的 math.deg 有问题...如果有人知道修复方法，请清理此代码！
与此同时，这里是 LUA 中 Haversine 的实现（与 Redis 一起使用！）
function calcDist(lat1, lon1, lat2, lon2)
lat1= lat1*0.0174532925
lat2= lat2*0.0174532925
lon1= lon1*0.0174532925
lon2= lon2*0.0174532925
dlon = lon2-lon1
dlat = lat2-lat1
a = math.pow(math.sin(dlat/2),2) + math.cos(lat1) * math.cos(lat2) * math.pow(math.sin(dlon/2),2)
c = 2 * math.asin(math.sqrt(a))
dist = 6371 * c -- multiply by 0.621371 to convert to miles
return dist
end
干杯！
【讨论】：
不要忘记将变量 (dlon,dlat,a,c,dist) 声明为local，以免污染全局空间。【参考方案7】：
FSharp 版本，使用里程：
let radialDistanceHaversine location1 location2 : float =
let degreeToRadian degrees = degrees * System.Math.PI / 180.0
let earthRadius = 3959.0
let deltaLat = location2.Latitude - location1.Latitude |> degreeToRadian
let deltaLong = location2.Longitude - location1.Longitude |> degreeToRadian
let a =
(deltaLat / 2.0 |> sin) ** 2.0
+ (location1.Latitude |> degreeToRadian |> cos)
* (location2.Latitude |> degreeToRadian |> cos)
* (deltaLong / 2.0 |> sin) ** 2.0
atan2 (a |> sqrt) (1.0 - a |> sqrt)
* 2.0
* earthRadius
【讨论】：
【参考方案8】：飞镖语言：
import 'dart:math' show cos, sqrt, asin;
double calculateDistance(LatLng l1, LatLng l2)
const p = 0.017453292519943295;
final a = 0.5 -
cos((l2.latitude - l1.latitude) * p) / 2 +
cos(l1.latitude * p) *
cos(l2.latitude * p) *
(1 - cos((l2.longitude - l1.longitude) * p)) /
2;
return 12742 * asin(sqrt(a));
【讨论】：
【参考方案9】：准确计算经纬度点之间的距离所需的函数很复杂，陷阱也很多。由于存在很大的不准确性（地球不是完美的球体），我不会推荐半正弦或其他球形解决方案。 vincenty formula 更好，但在某些情况下会抛出错误，即使编码正确。
我建议不要自己编写函数，而是使用 geopy，它已经实现了非常准确的 geographiclib 用于距离计算 (paper from author)。
#pip install geopy
from geopy.distance import geodesic
NY = [40.71278,-74.00594]
Beijing = [39.90421,116.40739]
print("WGS84: ",geodesic(NY, Beijing).km) #WGS84 is Standard
print("Intl24: ",geodesic(NY, Beijing, ellipsoid='Intl 1924').km) #geopy includes different ellipsoids
print("Custom ellipsoid: ",geodesic(NY, Beijing, ellipsoid=(6377., 6356., 1 / 297.)).km) #custom ellipsoid
#supported ellipsoids:
#model major (km) minor (km) flattening
#'WGS-84': (6378.137, 6356.7523142, 1 / 298.257223563)
#'GRS-80': (6378.137, 6356.7523141, 1 / 298.257222101)
#'Airy (1830)': (6377.563396, 6356.256909, 1 / 299.3249646)
#'Intl 1924': (6378.388, 6356.911946, 1 / 297.0)
#'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465)
#'GRS-67': (6378.1600, 6356.774719, 1 / 298.25)
这个库的唯一缺点是它不支持向量化计算。 对于矢量化计算，您可以使用新的geovectorslib。
#pip install geovectorslib
from geovectorslib import inverse
print(inverse(lats1,lons1,lats2,lons2)['s12'])
lats 和 lons 是 numpy 数组。 Geovectorslib 非常准确且速度极快！不过，我还没有找到改变椭球的解决方案。 WGS84 椭球作为标准使用，这是大多数用途的最佳选择。
【讨论】：
【参考方案10】：这是 Erlang 的实现
lat_lng(Lat1, Lon1=_Point1, Lat2, Lon2=_Point2) ->
P = math:pi() / 180,
R = 6371, % Radius of Earth in KM
A = 0.5 - math:cos((Lat2 - Lat1) * P) / 2 +
math:cos(Lat1 * P) * math:cos(Lat2 * P) * (1 - math:cos((Lon2 - Lon1) * P))/2,
R * 2 * math:asin(math:sqrt(A)).
【讨论】：
【参考方案11】：这是一个简单的 javascript 函数，可能对 link.. 有用，但我们使用的是 google earth javascript 插件而不是地图
function getApproximateDistanceUnits(point1, point2)
var xs = 0;
var ys = 0;
xs = point2.getX() - point1.getX();
xs = xs * xs;
ys = point2.getY() - point1.getY();
ys = ys * ys;
return Math.sqrt(xs + ys);
单位不是距离，而是相对于坐标的比率。还有其他相关的计算可以替代 getApproximateDistanceUnits 函数link here
然后我用这个函数来查看一个纬度经度是否在半径之内
function isMapPlacemarkInRadius(point1, point2, radi)
if (point1 && point2)
return getApproximateDistanceUnits(point1, point2) <= radi;
else
return 0;
点可以定义为
$$.getPoint = function(lati, longi)
var location =
x: 0,
y: 0,
getX: function() return location.x; ,
getY: function() return location.y;
;
location.x = lati;
location.y = longi;
return location;
;
然后你可以做你的事情来看看一个点是否在一个半径区域内：
//put it on the map if within the range of a specified radi assuming 100,000,000 units
var iconpoint = Map.getPoint(pp.latitude, pp.longitude);
var centerpoint = Map.getPoint(Settings.CenterLatitude, Settings.CenterLongitude);
//approx ~200 units to show only half of the globe from the default center radius
if (isMapPlacemarkInRadius(centerpoint, iconpoint, 120))
addPlacemark(pp.latitude, pp.longitude, pp.name);
else
otherSidePlacemarks.push(
latitude: pp.latitude,
longitude: pp.longitude,
name: pp.name
);
【讨论】：
【参考方案12】：//JAVA
public Double getDistanceBetweenTwoPoints(Double latitude1, Double longitude1, Double latitude2, Double longitude2)
final int RADIUS_EARTH = 6371;
double dLat = getRad(latitude2 - latitude1);
double dLong = getRad(longitude2 - longitude1);
double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(getRad(latitude1)) * Math.cos(getRad(latitude2)) * Math.sin(dLong / 2) * Math.sin(dLong / 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return (RADIUS_EARTH * c) * 1000;
private Double getRad(Double x)
return x * Math.PI / 180;
【讨论】：
【参考方案13】：我创建了这个小的 Javascript LatLng 对象，可能对某人有用。
var latLng1 = new LatLng(5, 3);
var latLng2 = new LatLng(6, 7);
var distance = latLng1.distanceTo(latLng2);
代码：
/**
* latLng point
* @param Number lat
* @param Number lng
* @returns LatLng
* @constructor
*/
function LatLng(lat,lng)
this.lat = parseFloat(lat);
this.lng = parseFloat(lng);
this.__cache = ;
LatLng.prototype =
toString: function()
return [this.lat, this.lng].join(",");
,
/**
* calculate distance in km to another latLng, with caching
* @param LatLng latLng
* @returns Number distance in km
*/
distanceTo: function(latLng)
var cacheKey = latLng.toString();
if(cacheKey in this.__cache)
return this.__cache[cacheKey];
// the fastest way to calculate the distance, according to this jsperf test;
// http://jsperf.com/haversine-salvador/8
// http://***.com/questions/27928
var deg2rad = 0.017453292519943295; // === Math.PI / 180
var lat1 = this.lat * deg2rad;
var lng1 = this.lng * deg2rad;
var lat2 = latLng.lat * deg2rad;
var lng2 = latLng.lng * deg2rad;
var a = (
(1 - Math.cos(lat2 - lat1)) +
(1 - Math.cos(lng2 - lng1)) * Math.cos(lat1) * Math.cos(lat2)
) / 2;
var distance = 12742 * Math.asin(Math.sqrt(a)); // Diameter of the earth in km (2 * 6371)
// cache the distance
this.__cache[cacheKey] = distance;
return distance;
;
【讨论】：
【参考方案14】：您可以使用Haversine公式计算它：
a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
c = 2 ⋅ atan2( √a, √(1−a) )
d = R ⋅ c
下面给出了一个计算两点之间距离的例子
假设我必须计算新德里到伦敦之间的距离，那么我该如何使用这个公式：
New delhi co-ordinates= 28.7041° N, 77.1025° E
London co-ordinates= 51.5074° N, 0.1278° W
var R = 6371e3; // metres
var φ1 = 28.7041.toRadians();
var φ2 = 51.5074.toRadians();
var Δφ = (51.5074-28.7041).toRadians();
var Δλ = (0.1278-77.1025).toRadians();
var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
Math.cos(φ1) * Math.cos(φ2) *
Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c; // metres
d = d/1000; // km
【讨论】：
【参考方案15】：如果你使用的是python； 点安装geopy
from geopy.distance import geodesic
origin = (30.172705, 31.526725) # (latitude, longitude) don't confuse
destination = (30.288281, 31.732326)
print(geodesic(origin, destination).meters) # 23576.805481751613
print(geodesic(origin, destination).kilometers) # 23.576805481751613
print(geodesic(origin, destination).miles) # 14.64994773134371
【讨论】：
【参考方案16】：function distance($lat1, $lon1, $lat2, $lon2)
$pi80 = M_PI / 180;
$lat1 *= $pi80; $lon1 *= $pi80; $lat2 *= $pi80; $lon2 *= $pi80;
$dlat = $lat2 - $lat1;
$dlon = $lon2 - $lon1;
$a = sin($dlat / 2) * sin($dlat / 2) + cos($lat1) * cos($lat2) * sin($dlon / 2) * sin($dlon / 2);
$km = 6372.797 * 2 * atan2(sqrt($a), sqrt(1 - $a));
return $km;
【讨论】：
【参考方案17】：您也可以使用 geolib 之类的模块：
如何安装：
$ npm install geolib
使用方法：
import getDistance from 'geolib'
const distance = getDistance(
latitude: 51.5103, longitude: 7.49347 ,
latitude: "51° 31' N", longitude: "7° 28' E"
)
console.log(distance)
文档： https://www.npmjs.com/package/geolib
【讨论】：
【参考方案18】：这是一个 Scala 实现：
def calculateHaversineDistance(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double =
val long2 = lon2 * math.Pi / 180
val lat2 = lat2 * math.Pi / 180
val long1 = lon1 * math.Pi / 180
val lat1 = lat1 * math.Pi / 180
val dlon = long2 - long1
val dlat = lat2 - lat1
val a = math.pow(math.sin(dlat / 2), 2) + math.cos(lat1) * math.cos(lat2) * math.pow(math.sin(dlon / 2), 2)
val c = 2 * math.atan2(Math.sqrt(a), math.sqrt(1 - a))
val haversineDistance = 3961 * c // 3961 = radius of earth in miles
haversineDistance
【讨论】：
以上是关于计算两个经纬度点之间的距离？ （Haversine 公式）的主要内容，如果未能解决你的问题，请参考以下文章