计算2个GPS坐标之间的距离

如何计算两个GPS坐标之间的距离(使用经纬度)?

516774 次浏览
小开

我猜你想让它沿着地球的曲率运动。你的两点和地心在一个平面上。地球的中心是这个平面上的圆心，这两个点(大致)在这个圆的周长上。由此你可以通过求一点到另一点的角度来计算距离。

如果点的高度不一样，或者如果你需要考虑地球不是一个完美的球体，这就有点困难了。

小开
最佳答案

通过经纬度计算两个坐标之间的距离，包括一个Javascript实现。

西和南位置为负。 记住，分和秒不是60度，所以S31 30'是-31.50度

别忘了将角度转换为弧度。许多语言都有这个功能。或者它是一个简单的计算:radians = degrees * PI / 180。

function degreesToRadians(degrees) {
return degrees * Math.PI / 180;
}


function distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2) {
var earthRadiusKm = 6371;


var dLat = degreesToRadians(lat2-lat1);
var dLon = degreesToRadians(lon2-lon1);


lat1 = degreesToRadians(lat1);
lat2 = degreesToRadians(lat2);


var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
return earthRadiusKm * c;
}

下面是一些用法的例子:

distanceInKmBetweenEarthCoordinates(0,0,0,0)  // Distance between same
// points should be 0
0


distanceInKmBetweenEarthCoordinates(51.5, 0, 38.8, -77.1) // From London
// to Arlington
5918.185064088764
小开

这个Lua代码改编自维基百科和Robert Lipe的GPSbabel工具:

local EARTH_RAD = 6378137.0
-- earth's radius in meters (official geoid datum, not 20,000km / pi)


local radmiles = EARTH_RAD*100.0/2.54/12.0/5280.0;
-- earth's radius in miles


local multipliers = {
radians = 1, miles = radmiles, mi = radmiles, feet = radmiles * 5280,
meters = EARTH_RAD, m = EARTH_RAD, km = EARTH_RAD / 1000,
degrees = 360 / (2 * math.pi), min = 60 * 360 / (2 * math.pi)
}


function gcdist(pt1, pt2, units) -- return distance in radians or given units
--- this formula works best for points close together or antipodal
--- rounding error strikes when distance is one-quarter Earth's circumference
--- (ref: wikipedia Great-circle distance)
if not pt1.radians then pt1 = rad(pt1) end
if not pt2.radians then pt2 = rad(pt2) end
local sdlat = sin((pt1.lat - pt2.lat) / 2.0);
local sdlon = sin((pt1.lon - pt2.lon) / 2.0);
local res = sqrt(sdlat * sdlat + cos(pt1.lat) * cos(pt2.lat) * sdlon * sdlon);
res = res > 1 and 1 or res < -1 and -1 or res
res = 2 * asin(res);
if units then return res * assert(multipliers[units])
else return res
end
end
小开

这取决于你需要它有多准确。如果你需要精确的精度，最好看看使用椭球而不是球体的算法，比如Vincenty的算法，它精确到毫米。

小开

在SQL Server 2008中使用地理类型非常容易做到这一点。

SELECT geography::Point(lat1, lon1, 4326).STDistance(geography::Point(lat2, lon2, 4326))
-- computes distance in meters using eliptical model, accurate to the mm

4326是WGS84椭球地球模型的SRID

小开

寻找带谷歌的哈弗辛;以下是我的解决方案:

#include <math.h>
#include "haversine.h"


#define d2r (M_PI / 180.0)


//calculate haversine distance for linear distance
double haversine_km(double lat1, double long1, double lat2, double long2)
{
double dlong = (long2 - long1) * d2r;
double dlat = (lat2 - lat1) * d2r;
double a = pow(sin(dlat/2.0), 2) + cos(lat1*d2r) * cos(lat2*d2r) * pow(sin(dlong/2.0), 2);
double c = 2 * atan2(sqrt(a), sqrt(1-a));
double d = 6367 * c;


return d;
}


double haversine_mi(double lat1, double long1, double lat2, double long2)
{
double dlong = (long2 - long1) * d2r;
double dlat = (lat2 - lat1) * d2r;
double a = pow(sin(dlat/2.0), 2) + cos(lat1*d2r) * cos(lat2*d2r) * pow(sin(dlong/2.0), 2);
double c = 2 * atan2(sqrt(a), sqrt(1-a));
double d = 3956 * c;


return d;
}
小开

下面是c#语言(用纬度和弧度表示):

double CalculateGreatCircleDistance(double lat1, double long1, double lat2, double long2, double radius)
{
return radius * Math.Acos(
Math.Sin(lat1) * Math.Sin(lat2)
+ Math.Cos(lat1) * Math.Cos(lat2) * Math.Cos(long2 - long1));
}

如果你的纬度和长度是用角度表示的，那么除以180/PI就可以转换成弧度。

小开

一个T-SQL函数，我用来根据中心的距离选择记录

Create Function  [dbo].[DistanceInMiles]
(  @fromLatitude float ,
@fromLongitude float ,
@toLatitude float,
@toLongitude float
)
returns float
AS
BEGIN
declare @distance float


select @distance = cast((3963 * ACOS(round(COS(RADIANS(90-@fromLatitude))*COS(RADIANS(90-@toLatitude))+
SIN(RADIANS(90-@fromLatitude))*SIN(RADIANS(90-@toLatitude))*COS(RADIANS(@fromLongitude-@toLongitude)),15))
)as float)
return  round(@distance,1)
END
小开 
    private double deg2rad(double deg)
{
return (deg * Math.PI / 180.0);
}


private double rad2deg(double rad)
{
return (rad / Math.PI * 180.0);
}


private double GetDistance(double lat1, double lon1, double lat2, double lon2)
{
//code for Distance in Kilo Meter
double theta = lon1 - lon2;
double dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta));
dist = Math.Abs(Math.Round(rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344 * 1000, 0));
return (dist);
}


private double GetDirection(double lat1, double lon1, double lat2, double lon2)
{
//code for Direction in Degrees
double dlat = deg2rad(lat1) - deg2rad(lat2);
double dlon = deg2rad(lon1) - deg2rad(lon2);
double y = Math.Sin(dlon) * Math.Cos(lat2);
double x = Math.Cos(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) - Math.Sin(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(dlon);
double direct = Math.Round(rad2deg(Math.Atan2(y, x)), 0);
if (direct < 0)
direct = direct + 360;
return (direct);
}


private double GetSpeed(double lat1, double lon1, double lat2, double lon2, DateTime CurTime, DateTime PrevTime)
{
//code for speed in Kilo Meter/Hour
TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
double TimeDifferenceInSeconds = Math.Round(TimeDifference.TotalSeconds, 0);
double theta = lon1 - lon2;
double dist = Math.Sin(deg2rad(lat1)) * Math.Sin(deg2rad(lat2)) + Math.Cos(deg2rad(lat1)) * Math.Cos(deg2rad(lat2)) * Math.Cos(deg2rad(theta));
dist = rad2deg(Math.Acos(dist)) * 60 * 1.1515 * 1.609344;
double Speed = Math.Abs(Math.Round((dist / Math.Abs(TimeDifferenceInSeconds)) * 60 * 60, 0));
return (Speed);
}


private double GetDuration(DateTime CurTime, DateTime PrevTime)
{
//code for speed in Kilo Meter/Hour
TimeSpan TimeDifference = CurTime.Subtract(PrevTime);
double TimeDifferenceInSeconds = Math.Abs(Math.Round(TimeDifference.TotalSeconds, 0));
return (TimeDifferenceInSeconds);
}
小开

如果你需要更精确的值，可以使用看看这个。

文森特提公式是大地测量学中使用的两个相关的迭代方法 计算a曲面上两点之间的距离 球形，由Thaddeus Vincenty (1975a)开发 假设地球的形状是一个扁球体，以及 因此比大圆距离等方法更准确

第一个(直接)方法计算a点的位置 给定到另一点的距离和方位(方向)。第二个 (逆)法计算地理距离和方位角 在两个给定点之间。它们在大地测量中得到了广泛应用 因为它们精确到地球上0.5毫米(0.020″) 椭球。< / p >
小开

你可以在fssnip的f#中找到它的实现(有一些很好的解释)

以下是重要的部分:


let GreatCircleDistance<[<Measure>] 'u> (R : float<'u>) (p1 : Location) (p2 : Location) =
let degToRad (x : float<deg>) = System.Math.PI * x / 180.0<deg/rad>


let sq x = x * x
// take the sin of the half and square the result
let sinSqHf (a : float<rad>) = (System.Math.Sin >> sq) (a / 2.0<rad>)
let cos (a : float<deg>) = System.Math.Cos (degToRad a / 1.0<rad>)


let dLat = (p2.Latitude - p1.Latitude) |> degToRad
let dLon = (p2.Longitude - p1.Longitude) |> degToRad


let a = sinSqHf dLat + cos p1.Latitude * cos p2.Latitude * sinSqHf dLon
let c = 2.0 * System.Math.Atan2(System.Math.Sqrt(a), System.Math.Sqrt(1.0-a))


R * c
小开

c#版本的Haversine

double _eQuatorialEarthRadius = 6378.1370D;
double _d2r = (Math.PI / 180D);


private int HaversineInM(double lat1, double long1, double lat2, double long2)
{
return (int)(1000D * HaversineInKM(lat1, long1, lat2, long2));
}


private double HaversineInKM(double lat1, double long1, double lat2, double long2)
{
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r) * Math.Pow(Math.Sin(dlong / 2D), 2D);
double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;


return d;
}

这是一个。net小提琴，所以你可以用你自己的Lat/ long测试它。

小开

这是“Henry Vilinskiy”为MySQL和km改编的版本:

CREATE FUNCTION `CalculateDistanceInKm`(
fromLatitude float,
fromLongitude float,
toLatitude float,
toLongitude float
) RETURNS float
BEGIN
declare distance float;


select
6367 * ACOS(
round(
COS(RADIANS(90-fromLatitude)) *
COS(RADIANS(90-toLatitude)) +
SIN(RADIANS(90-fromLatitude)) *
SIN(RADIANS(90-toLatitude)) *
COS(RADIANS(fromLongitude-toLongitude))
,15)
)
into distance;


return  round(distance,3);
END;
小开

我需要在PowerShell中实现这个，希望它可以帮助其他人。 关于这个方法的一些注意事项

  1. 不要分割任何一行，否则计算就会出错
  2. 要以KM为单位计算，在$distance的计算中去掉* 1000
  3. 更改$earthsRadius = 3963.19059，并在$distance的计算中删除* 1000，以英里为单位计算距离
  4. 我用的是Haversine，因为其他帖子指出Vincenty的公式更准确

    Function MetresDistanceBetweenTwoGPSCoordinates($latitude1, $longitude1, $latitude2, $longitude2)
{
$Rad = ([math]::PI / 180);


$earthsRadius = 6378.1370 # Earth's Radius in KM
$dLat = ($latitude2 - $latitude1) * $Rad
$dLon = ($longitude2 - $longitude1) * $Rad
$latitude1 = $latitude1 * $Rad
$latitude2 = $latitude2 * $Rad


$a = [math]::Sin($dLat / 2) * [math]::Sin($dLat / 2) + [math]::Sin($dLon / 2) * [math]::Sin($dLon / 2) * [math]::Cos($latitude1) * [math]::Cos($latitude2)
$c = 2 * [math]::ATan2([math]::Sqrt($a), [math]::Sqrt(1-$a))


$distance = [math]::Round($earthsRadius * $c * 1000, 0) #Multiple by 1000 to get metres


Return $distance
}
小开

基于Roman Makarov对这个线程的回复的Java版本的Haversine算法

public class HaversineAlgorithm {


static final double _eQuatorialEarthRadius = 6378.1370D;
static final double _d2r = (Math.PI / 180D);


public static int HaversineInM(double lat1, double long1, double lat2, double long2) {
return (int) (1000D * HaversineInKM(lat1, long1, lat2, long2));
}


public static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
* Math.pow(Math.sin(dlong / 2D), 2D);
double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;


return d;
}


}
小开

下面是我在Python中使用的Haversine函数:

from math import pi,sqrt,sin,cos,atan2


def haversine(pos1, pos2):
lat1 = float(pos1['lat'])
long1 = float(pos1['long'])
lat2 = float(pos2['lat'])
long2 = float(pos2['long'])


degree_to_rad = float(pi / 180.0)


d_lat = (lat2 - lat1) * degree_to_rad
d_long = (long2 - long1) * degree_to_rad


a = pow(sin(d_lat / 2), 2) + cos(lat1 * degree_to_rad) * cos(lat2 * degree_to_rad) * pow(sin(d_long / 2), 2)
c = 2 * atan2(sqrt(a), sqrt(1 - a))
km = 6367 * c
mi = 3956 * c


return {"km":km, "miles":mi}
小开

一、关于“面包屑”方法

  1. 地球半径在不同的纬度上是不同的。在Haversine算法中必须考虑到这一点。
  2. 考虑轴承的变化，它将直线变成拱门(更长的)
  3. 考虑到速度变化将把拱门变成螺旋(比拱门更长或更短)
  4. 高度变化将使平面螺旋变成3D螺旋(再次变长)。这对丘陵地区非常重要。

下面是考虑#1和#2的C语言函数:

double   calcDistanceByHaversine(double rLat1, double rLon1, double rHeading1,
double rLat2, double rLon2, double rHeading2){
double rDLatRad = 0.0;
double rDLonRad = 0.0;
double rLat1Rad = 0.0;
double rLat2Rad = 0.0;
double a = 0.0;
double c = 0.0;
double rResult = 0.0;
double rEarthRadius = 0.0;
double rDHeading = 0.0;
double rDHeadingRad = 0.0;


if ((rLat1 < -90.0) || (rLat1 > 90.0) || (rLat2 < -90.0) || (rLat2 > 90.0)
|| (rLon1 < -180.0) || (rLon1 > 180.0) || (rLon2 < -180.0)
|| (rLon2 > 180.0)) {
return -1;
};


rDLatRad = (rLat2 - rLat1) * DEGREE_TO_RADIANS;
rDLonRad = (rLon2 - rLon1) * DEGREE_TO_RADIANS;
rLat1Rad = rLat1 * DEGREE_TO_RADIANS;
rLat2Rad = rLat2 * DEGREE_TO_RADIANS;


a = sin(rDLatRad / 2) * sin(rDLatRad / 2) + sin(rDLonRad / 2) * sin(
rDLonRad / 2) * cos(rLat1Rad) * cos(rLat2Rad);


if (a == 0.0) {
return 0.0;
}


c = 2 * atan2(sqrt(a), sqrt(1 - a));
rEarthRadius = 6378.1370 - (21.3847 * 90.0 / ((fabs(rLat1) + fabs(rLat2))
/ 2.0));
rResult = rEarthRadius * c;


// Chord to Arc Correction based on Heading changes. Important for routes with many turns and U-turns


if ((rHeading1 >= 0.0) && (rHeading1 < 360.0) && (rHeading2 >= 0.0)
&& (rHeading2 < 360.0)) {
rDHeading = fabs(rHeading1 - rHeading2);
if (rDHeading > 180.0) {
rDHeading -= 180.0;
}
rDHeadingRad = rDHeading * DEGREE_TO_RADIANS;
if (rDHeading > 5.0) {
rResult = rResult * (rDHeadingRad / (2.0 * sin(rDHeadingRad / 2)));
} else {
rResult = rResult / cos(rDHeadingRad);
}
}
return rResult;
}

2有一种更简单的方法，效果很好。

按平均速度。

Trip_distance = Trip_average_speed * Trip_time

由于GPS速度是由多普勒效应检测的，与[Lon,Lat]没有直接关系，如果不是主要的距离计算方法，至少可以考虑作为次要的(备份或校正)。

小开

Scala版本

  def deg2rad(deg: Double) = deg * Math.PI / 180.0


def rad2deg(rad: Double) = rad / Math.PI * 180.0


def getDistanceMeters(lat1: Double, lon1: Double, lat2: Double, lon2: Double) = {
val theta = lon1 - lon2
val dist = Math.sin(deg2rad(lat1)) * Math.sin(deg2rad(lat2)) + Math.cos(deg2rad(lat1)) *
Math.cos(deg2rad(lat2)) * Math.cos(deg2rad(theta))
Math.abs(
Math.round(
rad2deg(Math.acos(dist)) * 60 * 1.1515 * 1.609344 * 1000)
)
}
小开

PHP版本:

(如果你的坐标已经是弧度，移除所有deg2rad()。)

$R = 6371; // km
$dLat = deg2rad($lat2-$lat1);
$dLon = deg2rad($lon2-$lon1);
$lat1 = deg2rad($lat1);
$lat2 = deg2rad($lat2);


$a = sin($dLat/2) * sin($dLat/2) +
sin($dLon/2) * sin($dLon/2) * cos($lat1) * cos($lat2);


$c = 2 * atan2(sqrt($a), sqrt(1-$a));
$d = $R * $c;
小开

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

function distance(lat1, lon1, lat2, lon2) {
var p = 0.017453292519943295;    // Math.PI / 180
var 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
def distance(lat1, lon1, lat2, lon2):
p = 0.017453292519943295
a = 0.5 - cos((lat2 - lat1) * p)/2 + cos(lat1 * p) * cos(lat2 * p) * (1 - cos((lon2 - lon1) * p)) / 2
return 12742 * asin(sqrt(a))

为了完整起见:半正矢 on wiki。

小开

如果你使用的是。net，不要重新启动轮子。看到System.Device.Location。在另一个答案的注释中归功于fnx。

using System.Device.Location;


double lat1 = 45.421527862548828D;
double long1 = -75.697189331054688D;
double lat2 = 53.64135D;
double long2 = -113.59273D;


GeoCoordinate geo1 = new GeoCoordinate(lat1, long1);
GeoCoordinate geo2 = new GeoCoordinate(lat2, long2);


double distance = geo1.GetDistanceTo(geo2);
小开

下面是答案中的Swift实现

func degreesToRadians(degrees: Double) -> Double {
return degrees * Double.pi / 180
}


func distanceInKmBetweenEarthCoordinates(lat1: Double, lon1: Double, lat2: Double, lon2: Double) -> Double {


let earthRadiusKm: Double = 6371


let dLat = degreesToRadians(degrees: lat2 - lat1)
let dLon = degreesToRadians(degrees: lon2 - lon1)


let lat1 = degreesToRadians(degrees: lat1)
let lat2 = degreesToRadians(degrees: lat2)


let a = sin(dLat/2) * sin(dLat/2) +
sin(dLon/2) * sin(dLon/2) * cos(lat1) * cos(lat2)
let c = 2 * atan2(sqrt(a), sqrt(1 - a))
return earthRadiusKm * c
}
小开

我把上面的答案用在Scala程序中

import java.lang.Math.{atan2, cos, sin, sqrt}


def latLonDistance(lat1: Double, lon1: Double)(lat2: Double, lon2: Double): Double = {
val earthRadiusKm = 6371
val dLat = (lat2 - lat1).toRadians
val dLon = (lon2 - lon1).toRadians
val latRad1 = lat1.toRadians
val latRad2 = lat2.toRadians


val a = sin(dLat / 2) * sin(dLat / 2) + sin(dLon / 2) * sin(dLon / 2) * cos(latRad1) * cos(latRad2)
val c = 2 * atan2(sqrt(a), sqrt(1 - a))
earthRadiusKm * c
}

我对函数进行了curcurated，以便能够轻松地产生具有两个位置固定之一的函数，并且只需要一对lat/lon来产生距离。

小开

这是我在Elixir中的实现

defmodule Geo do
@earth_radius_km 6371
@earth_radius_sm 3958.748
@earth_radius_nm 3440.065
@feet_per_sm 5280


@d2r :math.pi / 180


def deg_to_rad(deg), do: deg * @d2r


def great_circle_distance(p1, p2, :km), do: haversine(p1, p2) * @earth_radius_km
def great_circle_distance(p1, p2, :sm), do: haversine(p1, p2) * @earth_radius_sm
def great_circle_distance(p1, p2, :nm), do: haversine(p1, p2) * @earth_radius_nm
def great_circle_distance(p1, p2, :m), do: great_circle_distance(p1, p2, :km) * 1000
def great_circle_distance(p1, p2, :ft), do: great_circle_distance(p1, p2, :sm) * @feet_per_sm


@doc """
Calculate the [Haversine](https://en.wikipedia.org/wiki/Haversine_formula)
distance between two coordinates. Result is in radians. This result can be
multiplied by the sphere's radius in any unit to get the distance in that unit.
For example, multiple the result of this function by the Earth's radius in
kilometres and you get the distance between the two given points in kilometres.
"""
def haversine({lat1, lon1}, {lat2, lon2}) do
dlat = deg_to_rad(lat2 - lat1)
dlon = deg_to_rad(lon2 - lon1)


radlat1 = deg_to_rad(lat1)
radlat2 = deg_to_rad(lat2)


a = :math.pow(:math.sin(dlat / 2), 2) +
:math.pow(:math.sin(dlon / 2), 2) *
:math.cos(radlat1) * :math.cos(radlat2)


2 * :math.atan2(:math.sqrt(a), :math.sqrt(1 - a))
end
end
小开

飞镖版本

半正矢算法。

import 'dart:math';


class GeoUtils {


static double _degreesToRadians(degrees) {
return degrees * pi / 180;
}


static double distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2) {
var earthRadiusKm = 6371;


var dLat = _degreesToRadians(lat2-lat1);
var dLon = _degreesToRadians(lon2-lon1);


lat1 = _degreesToRadians(lat1);
lat2 = _degreesToRadians(lat2);


var a = sin(dLat/2) * sin(dLat/2) +
sin(dLon/2) * sin(dLon/2) * cos(lat1) * cos(lat2);
var c = 2 * atan2(sqrt(a), sqrt(1-a));
return earthRadiusKm * c;
}
}
小开

我认为R中的一个算法版本仍然缺失:

gpsdistance<-function(lat1,lon1,lat2,lon2){


# internal function to change deg to rad


degreesToRadians<- function (degrees) {
return (degrees * pi / 180)
}


R<-6371e3  #radius of Earth in meters


phi1<-degreesToRadians(lat1) # latitude 1
phi2<-degreesToRadians(lat2) # latitude 2
lambda1<-degreesToRadians(lon1) # longitude 1
lambda2<-degreesToRadians(lon2) # longitude 2


delta_phi<-phi1-phi2 # latitude-distance
delta_lambda<-lambda1-lambda2 # longitude-distance


a<-sin(delta_phi/2)*sin(delta_phi/2)+
cos(phi1)*cos(phi2)*sin(delta_lambda/2)*
sin(delta_lambda/2)


cc<-2*atan2(sqrt(a),sqrt(1-a))


distance<- R * cc


return(distance)  # in meters
}
小开

下面是Kotlin的一个变种:

import kotlin.math.*


class HaversineAlgorithm {


companion object {
private const val MEAN_EARTH_RADIUS = 6371.008
private const val D2R = Math.PI / 180.0
}


private fun haversineInKm(lat1: Double, lon1: Double, lat2: Double, lon2: Double): Double {
val lonDiff = (lon2 - lon1) * D2R
val latDiff = (lat2 - lat1) * D2R
val latSin = sin(latDiff / 2.0)
val lonSin = sin(lonDiff / 2.0)
val a = latSin * latSin + (cos(lat1 * D2R) * cos(lat2 * D2R) * lonSin * lonSin)
val c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a))
return MEAN_EARTH_RADIUS * c
}
}
小开

对于java

public static double degreesToRadians(double degrees) {
return degrees * Math.PI / 180;
}


public static double distanceInKmBetweenEarthCoordinates(Location location1, Location location2) {
double earthRadiusKm = 6371;


double dLat = degreesToRadians(location2.getLatitude()-location1.getLatitude());
double dLon = degreesToRadians(location2.getLongitude()-location1.getLongitude());


double lat1 = degreesToRadians(location1.getLatitude());
double lat2 = degreesToRadians(location2.getLatitude());


double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
return earthRadiusKm * c;
}
小开

对于任何寻找Delphi/Pascal版本的人:

function GreatCircleDistance(const Lat1, Long1, Lat2, Long2: Double): Double;
var
Lat1Rad, Long1Rad, Lat2Rad, Long2Rad: Double;
const
EARTH_RADIUS_KM = 6378;
begin
Lat1Rad  := DegToRad(Lat1);
Long1Rad := DegToRad(Long1);
Lat2Rad  := DegToRad(Lat2);
Long2Rad := DegToRad(Long2);
Result   := EARTH_RADIUS_KM * ArcCos(Cos(Lat1Rad) * Cos(Lat2Rad) * Cos(Long1Rad - Long2Rad) + Sin(Lat1Rad) * Sin(Lat2Rad));
end;

我对这个代码没有任何功劳，我最初是在一个公共论坛上发现Gary William发布的。

小开

在Python中，你可以使用geopy库使用WGS84椭球来计算测地线距离:

from geopy.distance import geodesic
newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)
print(geodesic(newport_ri, cleveland_oh).km)
小开

Unity版本c#

半正矢算法。

public float Distance(float lat1, float lon1, float lat2, float lon2)
{
var earthRadiusKm = 6371;


var dLat = (lat2 - lat1) * Mathf.Rad2Deg;
var dLon = (lon2 - lon1) * Mathf.Rad2Deg;


var a = Mathf.Sin(dLat / 2) * Mathf.Sin(dLat / 2) +
Mathf.Sin(dLon / 2) * Mathf.Sin(dLon / 2) *
Mathf.Cos(lat1 * Mathf.Rad2Deg) * Mathf.Cos(lat2 * Mathf.Rad2Deg);


var c = 2 * Mathf.Atan2(Mathf.Sqrt(a), Mathf.Sqrt(1 - a));
return earthRadiusKm * c;
}
小开

打印稿版本

export const degreeToRadian = (degree: number) => {
return degree * Math.PI / 180;
}


export const distanceBetweenEarthCoordinatesInKm = (lat1: number, lon1: number, lat2: number, lon2: number) => {
const earthRadiusInKm = 6371;


const dLat = degreeToRadian(lat2 - lat1);
const dLon = degreeToRadian(lon2 - lon1);


lat1 = degreeToRadian(lat1);
lat2 = degreeToRadian(lat2);


const a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.sin(dLon / 2) * Math.sin(dLon / 2) * Math.cos(lat1) * Math.cos(lat2);
const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));


return earthRadiusInKm * c;
}

提交答案