打印稿版本

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

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

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

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

别忘了把角度转换成弧度。许多语言都有这个功能。或者它是一个简单的计算:弧度=角度* 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

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

以下是重要的部分:

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

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

    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

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

一、关于“面包屑”方法

地球半径在不同的纬度上是不同的。在Haversine算法中必须考虑到这一点。 考虑轴承的变化，它将直线变成拱门(更长的) 考虑到速度变化将把拱门变成螺旋(比拱门更长或更短) 高度变化将使平面螺旋变成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]没有直接关系，如果不是主要的距离计算方法，至少可以考虑作为次要的(备份或校正)。