计算两个经纬度点之间的距离？（半正弦式）

此链接可能对您有所帮助，因为它详细介绍了如何使用半正弦公式来计算距离。

摘录：

这个脚本[在 Javascript 中] 计算两点之间的大圆距离 – 也就是说，在地球表面的最短距离 – 使用 “半正弦”公式。

function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {
  var R = 6371; // Radius of the earth in km
  var dLat = deg2rad(lat2-lat1);  // deg2rad below
  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 km
  return d;
}

function deg2rad(deg) {
  return deg * (Math.PI/180)
}

要计算球体上两点之间的距离，您需要进行大圆计算。

如果您需要将距离重新投影到平面，有许多 C/C++ 库可以帮助 MapTools 进行地图投影。为此，您将需要各种坐标系的投影字符串。

您可能还会发现 MapWindow 是可视化点的有用工具。此外，作为它的开源，它是如何使用 proj.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;
}

纬度和经度为十进制。我没有使用 min（） 进行 asin（） 调用，因为我使用的距离太小了，以至于它们不需要它。

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

额外更新：

如果您使用的是 iOS4 或更高版本，那么 Apple 会提供一些方法来执行此操作，因此可以通过以下方式实现相同的功能：

-(double)kilometresBetweenPlace1:(CLLocationCoordinate2D) place1 andPlace2:(CLLocationCoordinate2D) place2 {

    MKMapPoint  start, finish;


    start = MKMapPointForCoordinate(place1);
    finish = MKMapPointForCoordinate(place2);

    return MKMetersBetweenMapPoints(start, finish) / 1000;
}

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

在MySQL中列出表中指定点之间的距离（我们使用随机点 - 纬度：45.20327，经度：23.7806）小于50公里的所有点（表字段为coord_lat和coord_long）：

列出所有具有 DISTANCE<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 distanta 
FROM obiective 
WHERE coord_lat<>'' 
    AND coord_long<>'' 
HAVING distanta<50 
ORDER BY distanta desc

上面的例子在MySQL 5.0.95和5.5.16（Linux）中进行了测试。

这是一个简单的PHP函数，它将给出一个非常合理的近似值（在+/-1%的误差范围内）。

<?php
function 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;
}
?>

如前所述;地球不是一个球体。这就像马克·麦格维尔（Mark McGwire）决定练习的老棒球一样——它充满了凹痕和颠簸。更简单的计算（像这样）将其视为一个球体。

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

欲了解更多信息： 维基百科地理距离

这是半正弦公式的 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));
}

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

这是一个简单的 javascript 函数，可能在此链接中有用。.不知何故相关，但我们使用的是 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 函数链接的其他相关计算

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

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

        }

这里有一个很好的例子来计算距离，用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，此实现将根据您传递的枚举值为您提供以 KM 或英里为单位的结果。

Public Enum DistanceType
    Miles
    KiloMeters
End Enum

Public Structure Position
    Public Latitude As Double
    Public Longitude As Double
End 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 Math
earth_radius_mi = 3959
radians = 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。

对于大多数情况来说，haversine 绝对是一个很好的公式，其他答案已经包含它，所以我不打算占用空间。但重要的是要注意，无论使用什么公式（是的，不仅仅是一个）。因为精度范围很大，而且需要计算时间。公式的选择需要更多的思考，而不是一个简单的不费吹灰之力的答案。

这篇来自美国宇航局的人的帖子是我在讨论这些选项时发现的最好的帖子

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

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

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 km
  var dLat = deg2rad(lat2-lat1);  // deg2rad below
  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; 
  var miles = d / 1.609344; 

if ( units == 'km' ) {  
return d; 
 } else {
return miles;
}}

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

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

function distance(lat1, lon1, lat2, lon2) {
  const r = 6371; // km
  const p = Math.PI / 180;

  const a = 0.5 - Math.cos((lat2 - lat1) * p) / 2
                + Math.cos(lat1 * p) * Math.cos(lat2 * p) *
                  (1 - Math.cos((lon2 - lon1) * p)) / 2;

  return 2 * r * Math.asin(Math.sqrt(a));
}

你可以玩我的 jsPerf 并在这里查看结果。

最近我需要在 python 中做同样的事情，所以这里有一个 python 实现：

from math import cos, asin, sqrt, pi

def distance(lat1, lon1, lat2, lon2):
    r = 6371 # km
    p = pi / 180

    a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p) * cos(lat2*p) * (1-cos((lon2-lon1)*p))/2
    return 2 * r * asin(sqrt(a))

为了完整起见：维基百科上的Haversine。

这是我在经过一些搜索后通过十进制度计算距离的 java 实现。我使用以公里为单位的世界平均半径（来自维基百科）。 İ如果您想要结果英里，请使用以英里为单位的世界半径。

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;
}

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

使用包中的函数计算两点之间的距离非常简单：distmgeosphere

distm(p1, p2, fun = distHaversine)

哪里：

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

由于地球不是完美的球形，椭球体的 Vincenty 公式可能是计算距离的最佳方法。因此，在您使用的包中：geosphere

distm(p1, p2, fun = distVincentyEllipsoid)

---

当然，你不一定要使用包，你也可以用一个函数计算基地的距离：geosphereR

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

在 Mysql 中使用以下函数，将参数传递为 usingPOINT(LONG,LAT)

CREATE FUNCTION `distance`(a POINT, b POINT)
 RETURNS double
    DETERMINISTIC
BEGIN

RETURN

GLength( LineString(( PointFromWKB(a)), (PointFromWKB(b)))) * 100000; -- To Make the distance in meters

END;

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

我不喜欢添加另一个答案，但 Google maps 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/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;
    }
};

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

干杯！

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

pip install haversine

Python 实现

原点是美国本土的中心。

from haversine import haversine, Unit
origin = (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);
}

这是半正弦公式的打字稿实现

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 km
    var 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 km
    return d;
}

可能有一个更简单、更正确的解决方案：地球的周长在赤道上是40,000公里，在格林威治（或任何经度）周期上大约是37,000公里。因此：

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

我同意它应该进行微调，因为我自己说过它是一个椭球体，所以要乘以余弦的半径会有所不同。但它更准确一些。与谷歌地图相比，它确实大大减少了错误。

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

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

a= 6378.137#equitorial radius in km
b= 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 lat1
    x1=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 lat2
    x2=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$
DECLARE
    v_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 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

如前所述，准确的计算应该考虑到地球不是一个完美的球体。以下是此处提供的各种算法的一些比较：

geoDistance(50,5,58,3)
Haversine: 899 km
Maymenn: 833 km
Keerthana: 897 km
google.maps.geometry.spherical.computeDistanceBetween(): 900 km

geoDistance(50,5,-58,-3)
Haversine: 12030 km
Maymenn: 11135 km
Keerthana: 10310 km
google.maps.geometry.spherical.computeDistanceBetween(): 12044 km

geoDistance(.05,.005,.058,.003)
Haversine: 0.9169 km
Maymenn: 0.851723 km
Keerthana: 0.917964 km
google.maps.geometry.spherical.computeDistanceBetween(): 0.917964 km

geoDistance(.05,80,.058,80.3)
Haversine: 33.37 km
Maymenn: 33.34 km
Keerthana: 33.40767 km
google.maps.geometry.spherical.computeDistanceBetween(): 33.40770 km

在很短的距离上，Keerthana的算法似乎与谷歌地图的算法一致。谷歌地图似乎没有遵循任何简单的算法，这表明它可能是这里最准确的方法。

无论如何，这是 Keerthana 算法的 Javascript 实现：

function geoDistance(lat1, lng1, lat2, lng2){
    const a = 6378.137; // equitorial radius in km
    const 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.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

这是以公里为单位计算距离的 SQL 实现，

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 HAVING
distance < 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 radians
    latPoint1 = 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 - 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);
    }
}

飞镖郎：

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 points
    on 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-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.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.47m
Mean absolute error: -2499.84m
Max relative error: 0.55%
Mean relative error: -0.02%

因此，在 100,000 对随机坐标上平均偏差 2.5 公里，这在大多数情况下可能是好的。

下一个选项是 Vincent 公式，它精确到毫米，具体取决于收敛标准，也可以矢量化。它确实存在对跖点附近的收敛问题。您可以通过放宽收敛条件使其在这些点收敛，但准确率会下降到 0.25% 或更高。在对跖点之外，Vincenty 将在平均小于 1.e-6 的相对误差内提供接近 Geographiclib 的结果。

这里提到的 Geographiclib 确实是当前的黄金标准。它有多种实现，而且速度相当快，特别是如果您使用的是 C++ 版本。

现在，如果您打算将 Python 用于 10k 点以上的任何内容，我建议您考虑我的矢量化实现。我根据自己的需要创建了一个带有矢量化 Vincenty 例程的 geovectorslib 库，该库使用 Geographiclib 作为近对足点的后备。以下是 100k 点与 Geographiclib 的比较。如您所见，它为逆法提供了高达 20 倍的改进，为 100k 点的直接方法提供了 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)

您可以使用半正弦公式计算它，即：

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

对于那些正在寻找基于 WGS-84 和 GRS-80 标准的 Excel 公式的人：

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

源

精确计算纬度点之间距离所需的函数很复杂，而且存在很多缺陷。我不推荐 haversine 或其他球形解决方案，因为有很大的不准确之处（地球不是一个完美的球体）。vincenty 公式更好，但在某些情况下会抛出错误，即使编码正确也是如此。

我建议使用 geopy，而不是自己编写函数，它已经实现了非常准确的 geographiclib 用于距离计算（作者的论文）。

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

lat 和 lons 是 numpy 数组。Geovectorslib 非常准确且速度极快！不过，我还没有找到更改椭球体的解决方案。WGS84椭球体作为标准使用，是大多数用途的最佳选择。

我在 R 中制作了一个自定义函数，使用 R 基础包中提供的函数计算两个空间点之间的半弦距离 （km）。

custom_hav_dist <- function(lat1, lon1, lat2, lon2) {
R <- 6371
Radian_factor <- 0.0174533
lat_1 <- (90-lat1)*Radian_factor
lat_2 <- (90-lat2)*Radian_factor
diff_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：要以英里为单位计算距离，请将函数 （6371） 中的 R 替换为 3958.756（对于海里，请使用 3440.065）。

如果您使用的是 python; pip 安装 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

如果您想要行驶距离/路线（在此处发布，因为这是谷歌上两点之间距离的第一个结果，但对于大多数人来说，行驶距离更有用），您可以使用谷歌地图距离矩阵服务：

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

您也可以使用像 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

下面是一个 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
  }

要计算由纬度和经度指定的两点之间的距离，可以使用半正弦公式或文森特公式。这两个公式都考虑了地球的曲率，并提供了相当准确的结果。以下是每种方法的简要说明：

半正弦公式： 半正弦公式计算假设一个完美球体（在本例中为地球）上两点之间的大圆距离。它比 Vincent 的公式更简单、更快速，但精度可能略低，尤其是对于更长的距离。

半正弦公式如下：

a = sin²(Δlat/2) + cos(lat1) * cos(lat2) * sin²(Δlon/2)
c = 2 * atan2(√a, √(1-a))
d = R * c

哪里：

Δlat is the difference in latitude between the two points.
Δlon is the difference in longitude between the two points.
lat1 and lat2 are the latitudes of the two points.
R is the radius of the Earth (mean radius: 6,371 km).

您可以在半正弦实现中找到直接实现

文森特公式： Vincenty 的公式更复杂，但为计算地球椭球模型上两点之间的距离提供了更高的精度。它们考虑了地球的形状，并为地球上的任何一对点提供准确的结果。

Vincenty 公式有两种变体：直接公式，用于计算点之间的距离，以及逆公式，用于计算初始方位角、最终方位角和点之间的距离。

要使用 Vincent 的公式，您需要使用所选的编程语言实现特定的方程式。

直接配方： 直接公式在给定起点、初始方位角和距离的情况下计算目的地点（纬度和经度）。当您知道起点、所需距离和到达目的地的初始方位时，它很有用。

计算公式如下：

α1 = atan2((1 - f) * tan(lat1), cos(α0))
sinσ1 = (cos(U2) * sin(α1))^2 + (cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(α1))^2
sinσ = sqrt(sinσ1)
cosσ = sin(U1) * sin(U2) + cos(U1) * cos(U2) * cos(α1)
σ = atan2(sinσ, cosσ)
sinα = cos(U1) * cos(U2) * sin(α1) / sinσ
cos2αm = 1 - sinα^2
C = (f / 16) * cos2αm * (4 + f * (4 - 3 * cos2αm))
λ = L + (1 - C) * f * sinα * (σ + C * sinσ * (cos(2 * σ) + C * cosσ * (-1 + 2 * cos(2 * σ)^2)))
Δ = L2 - λ
L = λ

哪里：

lat1, lon1: Starting point latitude and longitude in radians.
α0: Initial bearing in radians.
U1 = atan((1 - f) * tan(lat1))
U2 = atan((1 - f) * tan(lat2)), where lat2 is calculated iteratively using the formula above.
f = (a - b) / a, where a and b are the equatorial and polar radii of the Earth, respectively.

逆公式： 逆公式计算地球表面两点之间的距离、初始方位角和最终方位角。当您知道两点的纬度和经度时，它很有用。

计算公式如下：

L = lon2 - lon1
λ = L
λʹ = 2π + atan2(y, x)
σ = atan2(yʹ, xʹ)
C = (f / 16) * cos^2αm * (4 + f * (4 - 3 * cos^2αm))
λ = L + (1 - C) * f * sinα * (σ + C * sinσ * (cos(2 * σ) + C * cosσ * (-1 + 2 * cos(2 * σ)^2)))

哪里：

lat1, lon1: Latitude and longitude of the first point in radians.
lat2, lon2: Latitude and longitude of the second point in radians.
L = lon2 - lon1
U1 = atan((1 - f) * tan(lat1))
U2 = atan((1 - f) * tan(lat2))
sinα = cos(U2) * sin(L)
cosα = sqrt(1 - sinα^2)
cos^2αm = cosα^2 * cosσ^2
x = σ - sinα * sinσ
y = λʹ - sinα * sinσ
xʹ = cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(L)
yʹ = cos(U2) * sin(L)