计算两个经纬度点之间的距离?(半正矢公式)

如何计算由经纬度指定的两点之间的距离?

为了澄清，我想用千米来表示距离;这些点使用WGS84系统，我想了解可用方法的相对准确性。

1082130 次浏览
小开
最佳答案

这链接可能对你有帮助，因为它详细说明了半正矢公式的使用来计算距离。

摘录:

这个脚本计算两点之间的大圆距离也就是说，在地球表面上的最短距离-使用“半正矢”公式。< / p >

function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {var R = 6371; // Radius of the earth in kmvar dLat = deg2rad(lat2-lat1);  // deg2rad belowvar dLon = deg2rad(lon2-lon1);var a =Math.sin(dLat/2) * Math.sin(dLat/2) +Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2)) *Math.sin(dLon/2) * Math.sin(dLon/2);var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));var d = R * c; // Distance in kmreturn d;}
function deg2rad(deg) {return deg * (Math.PI/180)}
小开

要计算球面上两点之间的距离，你需要做大圆计算。

如果你需要将距离重新投影到一个平面上，有许多C/ c++库可以帮助你在MapTools处进行地图投影。要做到这一点，你需要不同坐标系的投影字符串。

你可能还会发现MapWindow是一个可视化要点的有用工具。此外，由于它是开源的，它是如何使用project.dll库的有用指南，它似乎是核心的开源投影库。

小开

下面是一个c#实现:

static class DistanceAlgorithm{const double PIx = 3.141592653589793;const double RADIUS = 6378.16;
/// <summary>/// Convert degrees to Radians/// </summary>/// <param name="x">Degrees</param>/// <returns>The equivalent in radians</returns>public static double Radians(double x){return x * PIx / 180;}
/// <summary>/// Calculate the distance between two places./// </summary>/// <param name="lon1"></param>/// <param name="lat1"></param>/// <param name="lon2"></param>/// <param name="lat2"></param>/// <returns></returns>public static double DistanceBetweenPlaces(double lon1,double lat1,double lon2,double lat2){double dlon = Radians(lon2 - lon1);double dlat = Radians(lat2 - lat1);
double a = (Math.Sin(dlat / 2) * Math.Sin(dlat / 2)) + Math.Cos(Radians(lat1)) * Math.Cos(Radians(lat2)) * (Math.Sin(dlon / 2) * Math.Sin(dlon / 2));double angle = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));return angle * RADIUS;}
}
小开

非常感谢这一切。我在Objective-C iPhone应用程序中使用了以下代码:

const double PIx = 3.141592653589793;const double RADIO = 6371; // Mean radius of Earth in Km
double convertToRadians(double val) {
return val * PIx / 180;}
-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {
double dlon = convertToRadians(place2.longitude - place1.longitude);double dlat = convertToRadians(place2.latitude - place1.latitude);
double a = ( pow(sin(dlat / 2), 2) + cos(convertToRadians(place1.latitude))) * cos(convertToRadians(place2.latitude)) * pow(sin(dlon / 2), 2);double angle = 2 * asin(sqrt(a));
return angle * RADIO;}

纬度和经度是十进制的。我没有在asin()调用中使用min()，因为我使用的距离非常小，以至于它们不需要min()。

它给出了错误的答案，直到我传入弧度的值-现在它几乎与从苹果地图应用程序中获得的值相同:-)

额外的更新:

如果你使用的是iOS4或更高版本，那么苹果会提供一些方法来实现相同的功能:

-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {
MKMapPoint  start, finish;

start = MKMapPointForCoordinate(place1);finish = MKMapPointForCoordinate(place2);
return MKMetersBetweenMapPoints(start, finish) / 1000;}
小开

我在这里发布了我的工作示例。

列出表中指定点(我们使用一个随机点-纬度:45.20327，长:23.7806)之间距离小于50 KM的所有点，纬度和amp;经度，在MySQL(表字段是coord_lat和coord_long):

列出所有距离为50，单位为千米(考虑地球半径为6371千米):

SELECT denumire, (6371 * acos( cos( radians(45.20327) ) * cos( radians( coord_lat ) ) * cos( radians( 23.7806 ) - radians(coord_long) ) + sin( radians(45.20327) ) * sin( radians(coord_lat) ) )) AS distantaFROM obiectiveWHERE coord_lat<>''AND coord_long<>''HAVING distanta<50ORDER BY distanta desc

上面的例子是在MySQL 5.0.95和5.5.16 (Linux)中测试的。

小开

这是一个简单的PHP函数，它将给出一个非常合理的近似值(误差小于+/-1%)。

<?phpfunction distance($lat1, $lon1, $lat2, $lon2) {
$pi80 = M_PI / 180;$lat1 *= $pi80;$lon1 *= $pi80;$lat2 *= $pi80;$lon2 *= $pi80;
$r = 6372.797; // mean radius of Earth in km$dlat = $lat2 - $lat1;$dlon = $lon2 - $lon1;$a = sin($dlat / 2) * sin($dlat / 2) + cos($lat1) * cos($lat2) * sin($dlon / 2) * sin($dlon / 2);$c = 2 * atan2(sqrt($a), sqrt(1 - $a));$km = $r * $c;
//echo '<br/>'.$km;return $km;}?>

如前所述;地球不是一个球体。它就像马克·麦奎尔决定用来练习的一个很旧很旧的棒球——到处都是凹痕和凸起。简单的计算(像这样)把它当作一个球体。

不同的方法或多或少的精确取决于你在这个不规则的卵形上的位置以及你的点之间的距离(它们越近，绝对误差范围就越小)。你的期望越精确，计算就越复杂。

更多信息:地理距离

小开

下面是Haversine公式的java实现。

public final static double AVERAGE_RADIUS_OF_EARTH_KM = 6371;public int calculateDistanceInKilometer(double userLat, double userLng,double venueLat, double venueLng) {
double latDistance = Math.toRadians(userLat - venueLat);double lngDistance = Math.toRadians(userLng - venueLng);
double a = Math.sin(latDistance / 2) * Math.sin(latDistance / 2)+ Math.cos(Math.toRadians(userLat)) * Math.cos(Math.toRadians(venueLat))* Math.sin(lngDistance / 2) * Math.sin(lngDistance / 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return (int) (Math.round(AVERAGE_RADIUS_OF_EARTH_KM * c));}

请注意，这里我们将答案四舍五入到最近的km。

小开

这是一个简单的javascript函数，可能从链接有用。不知何故相关，但我们使用谷歌地球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函数链接在这里

然后我使用这个函数来查看经纬度是否在半径内

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 unitsvar 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 radiusif (isMapPlacemarkInRadius(centerpoint, iconpoint, 120)) {addPlacemark(pp.latitude, pp.longitude, pp.name);}else {otherSidePlacemarks.push({latitude: pp.latitude,longitude: pp.longitude,name: pp.name});
}
小开

这里有一个用PHP http://www.geodatasource.com/developers/php计算距离的好例子:

 function distance($lat1, $lon1, $lat2, $lon2, $unit) {
$theta = $lon1 - $lon2;$dist = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) +  cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta));$dist = acos($dist);$dist = rad2deg($dist);$miles = $dist * 60 * 1.1515;$unit = strtoupper($unit);
if ($unit == "K") {return ($miles * 1.609344);} else if ($unit == "N") {return ($miles * 0.8684);} else {return $miles;}}
小开

下面是VB的实现。NET，这个实现将根据您传递的Enum值以KM或Miles为单位给您结果。

Public Enum DistanceTypeMilesKiloMetersEnd Enum
Public Structure PositionPublic Latitude As DoublePublic Longitude As DoubleEnd Structure
Public Class Haversine
Public Function Distance(Pos1 As Position,Pos2 As Position,DistType As DistanceType) As Double
Dim R As Double = If((DistType = DistanceType.Miles), 3960, 6371)
Dim dLat As Double = Me.toRadian(Pos2.Latitude - Pos1.Latitude)
Dim dLon As Double = Me.toRadian(Pos2.Longitude - Pos1.Longitude)
Dim a As Double = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Cos(Me.toRadian(Pos1.Latitude)) * Math.Cos(Me.toRadian(Pos2.Latitude)) * Math.Sin(dLon / 2) * Math.Sin(dLon / 2)
Dim c As Double = 2 * Math.Asin(Math.Min(1, Math.Sqrt(a)))
Dim result As Double = R * c
Return result
End Function
Private Function toRadian(val As Double) As Double
Return (Math.PI / 180) * val
End Function
End Class
小开

我通过简化公式来简化计算。

下面是Ruby版本:

include Mathearth_radius_mi = 3959radians = lambda { |deg| deg * PI / 180 }coord_radians = lambda { |c| { :lat => radians[c[:lat]], :lng => radians[c[:lng]] } }
# from/to = { :lat => (latitude_in_degrees), :lng => (longitude_in_degrees) }def haversine_distance(from, to)from, to = coord_radians[from], coord_radians[to]cosines_product = cos(to[:lat]) * cos(from[:lat]) * cos(from[:lng] - to[:lng])sines_product = sin(to[:lat]) * sin(from[:lat])return earth_radius_mi * acos(cosines_product + sines_product)end
小开

你可以使用CLLocationDistance中的构建来计算这个:

CLLocation *location1 = [[CLLocation alloc] initWithLatitude:latitude1 longitude:longitude1];CLLocation *location2 = [[CLLocation alloc] initWithLatitude:latitude2 longitude:longitude2];[self distanceInMetersFromLocation:location1 toLocation:location2]
- (int)distanceInMetersFromLocation:(CLLocation*)location1 toLocation:(CLLocation*)location2 {CLLocationDistance distanceInMeters = [location1 distanceFromLocation:location2];return distanceInMeters;}

在你的例子中，如果你想要公里，只要除以1000。

小开

哈弗辛公式在大多数情况下都是很好的公式，其他答案已经包含了它所以我就不占用空间了。但重要的是要注意，无论使用什么公式(是的，不仅仅是一个)。因为可能的精度范围很大，以及所需的计算时间。公式的选择需要更多的思考，而不是简单的无脑答案。

这个帖子来自nasa的一个人，是我在讨论这些选项时发现的最好的一个

http://www.cs.nyu.edu/visual/home/proj/tiger/gisfaq.html

例如，如果您只是在100英里半径内按距离对行进行排序。地平公式比哈弗辛公式快得多。

HalfPi = 1.5707963;R = 3956; /* the radius gives you the measurement unit*/
a = HalfPi - latoriginrad;b = HalfPi - latdestrad;u = a * a + b * b;v = - 2 * a * b * cos(longdestrad - longoriginrad);c = sqrt(abs(u + v));return R * c;

注意这里只有一个余弦和一个平方根。在哈弗辛公式中有9个。

小开 
function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2,units) {var R = 6371; // Radius of the earth in kmvar dLat = deg2rad(lat2-lat1);  // deg2rad belowvar dLon = deg2rad(lon2-lon1);var a =Math.sin(dLat/2) * Math.sin(dLat/2) +Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2)) *Math.sin(dLon/2) * Math.sin(dLon/2);var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));var d = R * c;var miles = d / 1.609344;
if ( units == 'km' ) {return d;} else {return miles;}}

查克的解决方案，也适用于英里。

小开

在我的项目中，我需要计算很多点之间的距离，所以我继续尝试优化我在这里找到的代码。平均而言，在不同的浏览器中，我的新实现运行速度快2倍比得到最多好评的答案。

function distance(lat1, lon1, lat2, lon2) {var p = 0.017453292519943295;    // Math.PI / 180var c = Math.cos;var a = 0.5 - c((lat2 - lat1) * p)/2 +c(lat1 * p) * c(lat2 * p) *(1 - c((lon2 - lon1) * p))/2;
return 12742 * Math.asin(Math.sqrt(a)); // 2 * R; R = 6371 km}

您可以使用我的jsPerf并看到结果在这里。

最近我需要在python中做同样的事情，所以这里是python实现:

from math import cos, asin, sqrt, pi
def distance(lat1, lon1, lat2, lon2):p = pi/180a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p) * cos(lat2*p) * (1-cos((lon2-lon1)*p))/2return 12742 * asin(sqrt(a)) #2*R*asin...

为了完整起见:在维基百科上排名0。

小开

这是我的java实现计算距离经过一些搜索。我用的是世界平均半径(来自维基百科)，单位是千米。İf你想要的结果英里，然后使用世界半径英里。

public static double distanceLatLong2(double lat1, double lng1, double lat2, double lng2){double earthRadius = 6371.0d; // KM: use mile here if you want mile result
double dLat = toRadian(lat2 - lat1);double dLng = toRadian(lng2 - lng1);
double a = Math.pow(Math.sin(dLat/2), 2)  +Math.cos(toRadian(lat1)) * Math.cos(toRadian(lat2)) *Math.pow(Math.sin(dLng/2), 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
return earthRadius * c; // returns result kilometers}
public static double toRadian(double degrees){return (degrees * Math.PI) / 180.0d;}
小开

在其他答案中，中的实现缺失。

使用geosphere包中的distm函数计算两点之间的距离非常简单:

distm(p1, p2, fun = distHaversine)

地点:

p1 = longitude/latitude for point(s)p2 = longitude/latitude for point(s)# type of distance calculationfun = distCosine / distHaversine / distVincentySphere / distVincentyEllipsoid

由于地球不是完美的球形，所以第一个可能是计算距离的最佳方法。因此在你使用的geosphere包中，然后:

distm(p1, p2, fun = distVincentyEllipsoid)

当然你不一定要使用geosphere包，你也可以用一个函数来计算以R为基数的距离:

hav.dist <- function(long1, lat1, long2, lat2) {R <- 6371diff.long <- (long2 - long1)diff.lat <- (lat2 - lat1)a <- sin(diff.lat/2)^2 + cos(lat1) * cos(lat2) * sin(diff.long/2)^2b <- 2 * asin(pmin(1, sqrt(a)))d = R * breturn(d)}
小开

在Mysql中，使用下面的函数传递参数，使用POINT(LONG,LAT)

CREATE FUNCTION `distance`(a POINT, b POINT)RETURNS doubleDETERMINISTICBEGIN
RETURN
GLength( LineString(( PointFromWKB(a)), (PointFromWKB(b)))) * 100000; -- To Make the distance in meters
END;
小开 
//JAVApublic 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;}
小开

我不喜欢添加另一个答案，但谷歌地图API v.3具有球形几何(以及更多)。在将你的WGS84转换为十进制度后，你可以这样做:

<script src="http://maps.google.com/maps/api/js?sensor=false&libraries=geometry" type="text/javascript"></script>
distance = google.maps.geometry.spherical.computeDistanceBetween(new google.maps.LatLng(fromLat, fromLng),new google.maps.LatLng(toLat, toLng));

关于谷歌的计算有多精确，甚至使用了什么模型都没有任何消息(尽管它说的是“球面”而不是“大地水准面”。顺便说一下，“直线”距离显然不同于一个人在地球表面旅行的距离，而这似乎是每个人都在假设的。

小开

我已经创建了这个小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://stackoverflow.com/questions/27928var deg2rad = 0.017453292519943295; // === Math.PI / 180var 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 distancethis.__cache[cacheKey] = distance;
return distance;}};
小开

数学有问题，LUA的学位…如果有人知道修复，请清理这段代码!

与此同时，这里有一个Haversine在LUA中的实现(与Redis一起使用!)

function calcDist(lat1, lon1, lat2, lon2)lat1= lat1*0.0174532925lat2= lat2*0.0174532925lon1= lon1*0.0174532925lon2= lon2*0.0174532925
dlon = lon2-lon1dlat = 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 milesreturn distend

干杯!

小开 
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}));
小开

# 0

Python实现

原产地是美国毗连的中心。

from haversine import haversine, Unitorigin = (39.50, 98.35)paris = (48.8567, 2.3508)haversine(origin, paris, unit=Unit.MILES)

要得到以公里为单位的答案，只需设置unit=Unit.KILOMETERS(这是默认值)。

小开

下面是移植到Java的已接受的答案实现，以备任何人需要。

package com.project529.garage.util;

/*** Mean radius.*/private static double EARTH_RADIUS = 6371;
/*** Returns the distance between two sets of latitudes and longitudes in meters.* <p/>* Based from the following JavaScript SO answer:* http://stackoverflow.com/questions/27928/calculate-distance-between-two-latitude-longitude-points-haversine-formula,* which is based on https://en.wikipedia.org/wiki/Haversine_formula (error rate: ~0.55%).*/public double getDistanceBetween(double lat1, double lon1, double lat2, double lon2) {double dLat = toRadians(lat2 - lat1);double dLon = toRadians(lon2 - lon1);
double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) +Math.cos(toRadians(lat1)) * Math.cos(toRadians(lat2)) *Math.sin(dLon / 2) * Math.sin(dLon / 2);double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));double d = EARTH_RADIUS * c;
return d;}
public double toRadians(double degrees) {return degrees * (Math.PI / 180);}
小开

下面是哈弗辛公式的0号实现

static getDistanceFromLatLonInKm(lat1: number, lon1: number, lat2: number, lon2: number): number {var deg2Rad = deg => {return deg * Math.PI / 180;}
var r = 6371; // Radius of the earth in kmvar dLat = deg2Rad(lat2 - lat1);var dLon = deg2Rad(lon2 - lon1);var a =Math.sin(dLat / 2) * Math.sin(dLat / 2) +Math.cos(deg2Rad(lat1)) * Math.cos(deg2Rad(lat2)) *Math.sin(dLon / 2) * Math.sin(dLon / 2);var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));var d = r * c; // Distance in kmreturn d;}
小开

可能有一个更简单、更正确的解决方案:地球的周长在赤道上是40000公里，在格林威治(或任何经度)周期上约为37000公里。因此:

pythagoras = function (lat1, lon1, lat2, lon2) {function sqr(x) {return x * x;}function cosDeg(x) {return Math.cos(x * Math.PI / 180.0);}
var earthCyclePerimeter = 40000000.0 * cosDeg((lat1 + lat2) / 2.0);var dx = (lon1 - lon2) * earthCyclePerimeter / 360.0;var dy = 37000000.0 * (lat1 - lat2) / 360.0;
return Math.sqrt(sqr(dx) + sqr(dy));};

我同意它应该被微调，我自己说过它是一个椭球，所以半径乘以余弦值是不同的。但它更准确一点。与谷歌map相比，误差明显减小。

小开

在提供的代码中有一些错误，我在下面修复了它。

以上所有答案都假定地球是一个球体。然而，更精确的近似是扁球体。

a= 6378.137#equitorial radius in kmb= 6356.752#polar radius in km
def Distance(lat1, lons1, lat2, lons2):lat1=math.radians(lat1)lons1=math.radians(lons1)R1=(((((a**2)*math.cos(lat1))**2)+(((b**2)*math.sin(lat1))**2))/((a*math.cos(lat1))**2+(b*math.sin(lat1))**2))**0.5 #radius of earth at lat1x1=R1*math.cos(lat1)*math.cos(lons1)y1=R1*math.cos(lat1)*math.sin(lons1)z1=R1*math.sin(lat1)
lat2=math.radians(lat2)lons2=math.radians(lons2)R2=(((((a**2)*math.cos(lat2))**2)+(((b**2)*math.sin(lat2))**2))/((a*math.cos(lat2))**2+(b*math.sin(lat2))**2))**0.5 #radius of earth at lat2x2=R2*math.cos(lat2)*math.cos(lons2)y2=R2*math.cos(lat2)*math.sin(lons2)z2=R2*math.sin(lat2)    
return ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5
小开

下面是postgres SQL中的一个示例(以km为单位，用于英里版本，将1.609344替换为0.8684版本)

CREATE OR REPLACE FUNCTION public.geodistance(alat float, alng float, blat
float, blng  float)RETURNS float AS$BODY$DECLAREv_distance float;BEGIN
v_distance = asin( sqrt(sin(radians(blat-alat)/2)^2+ ((sin(radians(blng-alng)/2)^2) *cos(radians(alat)) *cos(radians(blat))))) * cast('7926.3352' as float) * cast('1.609344' as float) ;

RETURN v_distance;END$BODY$language plpgsql VOLATILE SECURITY DEFINER;alter function geodistance(alat float, alng float, blat float, blng float)owner to postgres;
小开

这个脚本[在PHP中]计算两点之间的距离。

public static function getDistanceOfTwoPoints($source, $dest, $unit='K') {$lat1 = $source[0];$lon1 = $source[1];$lat2 = $dest[0];$lon2 = $dest[1];
$theta = $lon1 - $lon2;$dist = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) +  cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta));$dist = acos($dist);$dist = rad2deg($dist);$miles = $dist * 60 * 1.1515;$unit = strtoupper($unit);
if ($unit == "K") {return ($miles * 1.609344);}else if ($unit == "M"){return ($miles * 1.609344 * 1000);}else if ($unit == "N") {return ($miles * 0.8684);}else {return $miles;}}
小开

下面是另一个转换为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 kilometerdLat = 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 kilometerend
小开

正如指出的那样，精确的计算应该考虑到地球不是一个完美的球体。以下是这里提供的各种算法的一些比较:

geoDistance(50,5,58,3)Haversine: 899 kmMaymenn: 833 kmKeerthana: 897 kmgoogle.maps.geometry.spherical.computeDistanceBetween(): 900 km
geoDistance(50,5,-58,-3)Haversine: 12030 kmMaymenn: 11135 kmKeerthana: 10310 kmgoogle.maps.geometry.spherical.computeDistanceBetween(): 12044 km
geoDistance(.05,.005,.058,.003)Haversine: 0.9169 kmMaymenn: 0.851723 kmKeerthana: 0.917964 kmgoogle.maps.geometry.spherical.computeDistanceBetween(): 0.917964 km
geoDistance(.05,80,.058,80.3)Haversine: 33.37 kmMaymenn: 33.34 kmKeerthana: 33.40767 kmgoogle.maps.geometry.spherical.computeDistanceBetween(): 33.40770 km

在小范围内，Keerthana的算法似乎与谷歌Maps的算法一致。谷歌Maps似乎没有遵循任何简单的算法，这表明它可能是这里最准确的方法。

不管怎样，这里是Keerthana算法的Javascript实现:

function geoDistance(lat1, lng1, lat2, lng2){const a = 6378.137; // equitorial radius in kmconst b = 6356.752; // polar radius in km
var sq = x => (x*x);var sqr = x => Math.sqrt(x);var cos = x => Math.cos(x);var sin = x => Math.sin(x);var radius = lat => sqr((sq(a*a*cos(lat))+sq(b*b*sin(lat)))/(sq(a*cos(lat))+sq(b*sin(lat))));
lat1 = lat1 * Math.PI / 180;lng1 = lng1 * Math.PI / 180;lat2 = lat2 * Math.PI / 180;lng2 = lng2 * Math.PI / 180;
var R1 = radius(lat1);var x1 = R1*cos(lat1)*cos(lng1);var y1 = R1*cos(lat1)*sin(lng1);var z1 = R1*sin(lat1);
var R2 = radius(lat2);var x2 = R2*cos(lat2)*cos(lng2);var y2 = R2*cos(lat2)*sin(lng2);var z2 = R2*sin(lat2);
return sqr(sq(x1-x2)+sq(y1-y2)+sq(z1-z2));}
小开

FSharp版本，使用里程:

let radialDistanceHaversine location1 location2 : float =let degreeToRadian degrees = degrees * System.Math.PI / 180.0let earthRadius = 3959.0let deltaLat = location2.Latitude - location1.Latitude |> degreeToRadianlet deltaLong = location2.Longitude - location1.Longitude |> degreeToRadianlet a =(deltaLat / 2.0 |> sin) ** 2.0+ (location1.Latitude |> degreeToRadian |> cos)* (location2.Latitude |> degreeToRadian |> cos)* (deltaLong / 2.0 |> sin) ** 2.0atan2 (a |> sqrt) (1.0 - a |> sqrt)* 2.0* earthRadius
小开

下面是SQL实现，以km为单位计算距离，

SELECT UserId, ( 3959 * acos( cos( radians( your latitude here ) ) * cos( radians(latitude) ) *cos( radians(longitude) - radians( your longitude here ) ) + sin( radians( your latitude here ) ) *sin( radians(latitude) ) ) ) AS distance FROM user HAVINGdistance < 5  ORDER BY distance LIMIT 0 , 5;

要了解编程语言实现的更多细节，您可以浏览在这里给出的php脚本

小开

Java实现在半正矢公式

double calculateDistance(double latPoint1, double lngPoint1,double latPoint2, double lngPoint2) {if(latPoint1 == latPoint2 && lngPoint1 == lngPoint2) {return 0d;}
final double EARTH_RADIUS = 6371.0; //km value;
//converting to radianslatPoint1 = Math.toRadians(latPoint1);lngPoint1 = Math.toRadians(lngPoint1);latPoint2 = Math.toRadians(latPoint2);lngPoint2 = Math.toRadians(lngPoint2);
double distance = Math.pow(Math.sin((latPoint2 - latPoint1) / 2.0), 2)+ Math.cos(latPoint1) * Math.cos(latPoint2)* Math.pow(Math.sin((lngPoint2 - lngPoint1) / 2.0), 2);distance = 2.0 * EARTH_RADIUS * Math.asin(Math.sqrt(distance));
return distance; //km value}
小开
计算距离——尤其是大距离——的主要挑战之一是解释地球的曲率。如果地球是平的，计算两点之间的距离就会像计算直线一样简单!哈弗辛公式包括一个常数(下面是R变量)，它表示地球的半径。根据你测量的是英里还是公里，它分别等于3956英里或6367公里基本公式是:

dlon = lon2 - lon1dlat = lat2 - lat1a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2c = 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 / 180var radlat2 = Math.PI * lat2 / 180var theta = lon1 - lon2var radtheta = Math.PI * theta / 180var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);dist = Math.acos(dist)dist = dist * 180 / Math.PIdist = dist * 60 * 1.1515if (unit == "K") {dist = dist * 1.609344}if (unit == "N") {dist = dist * 0.8684}return dist} catch (err) {console.log(err);}}
小开

飞镖lang:

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));}
小开

由于这是关于这个话题最受欢迎的讨论，我将在这里补充我从2019年底到2020年初的经验。为了补充现有的答案-我的重点是找到一个准确和快速(即向量化)的解决方案。

让我们从这里最常用的答案——哈弗辛方法开始。向量化是很简单的，参见下面python中的例子:

def haversine(lat1, lon1, lat2, lon2):"""Calculate the great circle distance between two pointson the earth (specified in decimal degrees)
All args must be of equal length.Distances are in meters.    
Ref:https://stackoverflow.com/questions/29545704/fast-haversine-approximation-python-pandashttps://ipython.readthedocs.io/en/stable/interactive/magics.html"""Radius = 6.371e6lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1dlat = 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 degreesy = 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%。我的实验显示偏差较小。以下是10万个随机点与我的库的比较，应该精确到毫米级:

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.2f}m".format(np.max(r1['s12']-r2['s12'])))print("Mean absolute error: {:4.2f}m".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.47mMean absolute error: -2499.84mMax relative error: 0.55%Mean relative error: -0.02%

因此，在10万对随机坐标上，平均偏差为2.5km，这可能对大多数情况都是好的。

下一个选择是Vincenty公式，精确到毫米，这取决于收敛标准，也可以向量化。它确实有在对跖点附近收敛的问题。你可以通过放宽收敛标准使其收敛于这些点，但准确度会下降到0.25%甚至更多。在对映点之外，Vincenty将提供与地理库相近的结果，相对误差小于1。平均是E-6。

这里提到的Geographiclib实际上是当前的黄金标准。它有几个实现，而且相当快，特别是如果你使用的是c++版本。

现在，如果您计划将Python用于任何超过10k点的内容，我建议考虑我的向量化实现。我根据自己的需要创建了一个带有向量化Vincenty例程的geovectorslib库，它使用地理库作为近对点的备份。下面是与地理学lib的100k点的比较。正如你所看到的，它为100k分提供了最多反向提升20倍，直接提升100倍的方法，差距将随着分数的增加而增加。准确地说，它将在1以内。e-5 rtol的地理库。

Direct method for 100,000 points94.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 points1.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)
小开

你可以用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° ELondon co-ordinates= 51.5074° N, 0.1278° W
var R = 6371e3; // metresvar φ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; // metresd = d/1000; // km
小开

对于那些寻找基于WGS-84 &grs - 80标准:

=ACOS(COS(RADIANS(90-Lat1))*COS(RADIANS(90-Lat2))+SIN(RADIANS(90-Lat1))*SIN(RADIANS(90-Lat2))*COS(RADIANS(Long1-Long2)))*6371

Source .

小开

精确计算中长点之间距离所需的函数是复杂的，陷阱也很多。我不推荐哈弗辛或其他球形的解决方案，因为有很大的不准确性(地球不是一个完美的球体)。vincenty公式更好，但在某些情况下会抛出错误，即使编码正确。

与其自己编写函数，我建议使用geopy，它已经实现了非常精确的geographiclib距离计算(作者论文)。

#pip install geopyfrom geopy.distance import geodesicNY = [40.71278,-74.00594]Beijing = [39.90421,116.40739]print("WGS84: ",geodesic(NY, Beijing).km) #WGS84 is Standardprint("Intl24: ",geodesic(NY, Beijing, ellipsoid='Intl 1924').km) #geopy includes different ellipsoidsprint("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 geovectorslibfrom geovectorslib import inverseprint(inverse(lats1,lons1,lats2,lons2)['s12'])

lat和lon是numpy数组。Geovectorslib是非常准确和非常快!我还没有找到改变椭球的方法。标准采用WGS84椭球，是大多数用途的最佳选择。

小开

我在R中做了一个自定义函数，使用R基本包中可用的函数来计算两个空间点之间的距离(km)。

custom_hav_dist <- function(lat1, lon1, lat2, lon2) {R <- 6371Radian_factor <- 0.0174533lat_1 <- (90-lat1)*Radian_factorlat_2 <- (90-lat2)*Radian_factordiff_long <-(lon1-lon2)*Radian_factor
distance_in_km <- 6371*acos((cos(lat_1)*cos(lat_2))+(sin(lat_1)*sin(lat_2)*cos(diff_long)))rm(lat1, lon1, lat2, lon2)return(distance_in_km)}

样例输出

custom_hav_dist(50.31,19.08,54.14,19.39)[1] 426.3987

PS:要计算以英里为单位的距离，请将函数R(6371)替换为3958.756(海里使用3440.065)。

小开

如果你正在使用python;PIP install geopy

from geopy.distance import geodesic

origin = (30.172705, 31.526725)  # (latitude, longitude) don't confusedestination = (30.288281, 31.732326)
print(geodesic(origin, destination).meters)  # 23576.805481751613print(geodesic(origin, destination).kilometers)  # 23.576805481751613print(geodesic(origin, destination).miles)  # 14.64994773134371
小开

如果你想要驾驶距离/路线(张贴在这里，因为这是谷歌上两点之间距离的第一个结果，但对大多数人来说，驾驶距离更有用)，你可以使用谷歌映射距离矩阵服务:

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});}}});});}
小开 
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;}
小开

下面是Erlang实现

lat_lng({Lat1, Lon1}=_Point1, {Lat2, Lon2}=_Point2) ->P = math:pi() / 180,R = 6371, % Radius of Earth in KMA = 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)).
小开

你也可以使用像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)
< p >文档:# 0 < / p >
小开

下面是一个Scala实现:

  def calculateHaversineDistance(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double = {val long2 = lon2 * math.Pi / 180val lat2 = lat2 * math.Pi / 180val long1 = lon1 * math.Pi / 180val lat1 = lat1 * math.Pi / 180
val dlon = long2 - long1val dlat = lat2 - lat1val 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 mileshaversineDistance}

提交答案