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


当前回答

一、关于“面包屑”方法

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

其他回答

一、关于“面包屑”方法

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

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测试它。

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;

这段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

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