一、关于“面包屑”方法

地球半径在不同的纬度上是不同的。在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]没有直接关系，如果不是主要的距离计算方法，至少可以考虑作为次要的(备份或校正)。

计算两个坐标之间的纬度和经度的距离，包括一个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

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

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

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

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

Unity版本c#

Haversine Algorithm。

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