为什么boost :: geometry地理Vincenty距离在赤道附近不正确？

Question

我需要一个函数来高精度地计算一对WGS 84位置之间的距离，并且我打算使用Boost几何中的geographic函数。

该升压几何设计合理的状态：

有Andoyer方法，既快速又精确，还有Vincenty方法，又慢又精确。

但是，在boost::geometry::distance同时使用Andoyer和Vincenty策略测试功能时，我得到了以下结果：

WGS 84 values (metres)
    Semimajor axis:         6378137.000000
    Flattening:             0.003353
    Semiminor axis:         6356752.314245

    Semimajor distance:     20037508.342789
    Semiminor distance:     19970326.371123

Boost geometry near poles
Andoyer function:
    Semimajor distance:     20037508.151445
    Semiminor distance:     20003917.164970
Vincenty function:
    Semimajor distance:     **19970326.180419**
    Semiminor distance:     20003931.266635

Boost geometry at poles
Andoyer function:
    Semimajor distance:     0.000000
    Semiminor distance:     0.000000
Vincenty function:
    Semimajor distance:     **19970326.371122**
    Semiminor distance:     20003931.458623

在Vincenty沿半长轴（即在赤道附近）距离小于周围的短半轴的距离北极和南极之间的轴。那是不正确的。

Semiminor和Andoyer距离看起来很合理。除非这些点位于地球的相对侧，否则boost Andoyer函数将返回零！

问题出在：Vincenty算法，其boost geometry实现还是我的测试代码？

测试代码：

/// boost geometry WGS84 distance issue

// Note: M_PI is not part of the C or C++ standards, _USE_MATH_DEFINES enables it
#define _USE_MATH_DEFINES
#include <boost/geometry.hpp>
#include <cmath>
#include <iostream>
#include <ios>

// WGS 84 parameters from: Eurocontrol WGS 84 Implementation Manual
// Version 2.4 Chapter 3, page 14

/// The Semimajor axis measured in metres.
/// This is the radius at the equator.
constexpr double a = 6378137.0;

/// Flattening, a ratio.
/// This is the flattening of the ellipse at the poles
constexpr double f = 1.0/298.257223563;

/// The Semiminor axis measured in metres.
/// This is the radius at the poles.
/// Note: this is derived from the Semimajor axis and the flattening.
/// See WGS 84 Implementation Manual equation B-2, page 69.
constexpr double b = a * (1.0 - f);

int main(int /*argc*/, char ** /*argv*/)
{
  std::cout.setf(std::ios::fixed);

  std::cout << "WGS 84 values (metres)\n";
  std::cout << "\tSemimajor axis:\t\t"   << a << "\n";
  std::cout << "\tFlattening:\t\t"       << f << "\n";
  std::cout << "\tSemiminor axis:\t\t"   << b << "\n\n";

  std::cout << "\tSemimajor distance:\t" << M_PI * a << "\n";
  std::cout << "\tSemiminor distance:\t" << M_PI * b << "\n";
  std::cout << std::endl;

  // Min value for delta. 0.000000014 causes Andoyer to fail.
  const double DELTA(0.000000015);

  // For boost::geometry:
  typedef boost::geometry::cs::geographic<boost::geometry::radian> Wgs84Coords;
  typedef boost::geometry::model::point<double, 2, Wgs84Coords> GeographicPoint;
  // Note boost points are Long & Lat NOT Lat & Long
  GeographicPoint near_north_pole   (0.0,  M_PI_2 - DELTA);
  GeographicPoint near_south_pole   (0.0, -M_PI_2 + DELTA);

  GeographicPoint near_equator_east ( M_PI_2 - DELTA, 0.0);
  GeographicPoint near_equator_west (-M_PI_2 + DELTA, 0.0);

  // Note: the default boost geometry spheroid is WGS84
  // #include <boost/geometry/core/srs.hpp>
  typedef boost::geometry::srs::spheroid<double> SpheroidType;
  SpheroidType spheriod;

  //#include <boost/geometry/strategies/geographic/distance_andoyer.hpp>
  typedef boost::geometry::strategy::distance::andoyer<SpheroidType>
                                                               AndoyerStrategy;
  AndoyerStrategy andoyer(spheriod);

  std::cout << "Boost geometry near poles\n";
  std::cout << "Andoyer function:\n";
  double andoyer_major(boost::geometry::distance(near_equator_east, near_equator_west, andoyer));
  std::cout << "\tSemimajor distance:\t" << andoyer_major << "\n";
  double andoyer_minor(boost::geometry::distance(near_north_pole, near_south_pole, andoyer));
  std::cout << "\tSemiminor distance:\t" << andoyer_minor << "\n";

  //#include <boost/geometry/strategies/geographic/distance_vincenty.hpp>
  typedef boost::geometry::strategy::distance::vincenty<SpheroidType>
                                                               VincentyStrategy;
  VincentyStrategy vincenty(spheriod);

  std::cout << "Vincenty function:\n";
  double vincenty_major(boost::geometry::distance(near_equator_east, near_equator_west, vincenty));
  std::cout << "\tSemimajor distance:\t" << vincenty_major << "\n";
  double vincenty_minor(boost::geometry::distance(near_north_pole, near_south_pole, vincenty));
  std::cout << "\tSemiminor distance:\t" << vincenty_minor << "\n\n";

  // Note boost points are Long & Lat NOT Lat & Long
  GeographicPoint north_pole   (0.0,  M_PI_2);
  GeographicPoint south_pole   (0.0, -M_PI_2);

  GeographicPoint equator_east ( M_PI_2, 0.0);
  GeographicPoint equator_west (-M_PI_2, 0.0);

  std::cout << "Boost geometry at poles\n";
  std::cout << "Andoyer function:\n";
  andoyer_major = boost::geometry::distance(equator_east, equator_west, andoyer);
  std::cout << "\tSemimajor distance:\t" << andoyer_major << "\n";
  andoyer_minor = boost::geometry::distance(north_pole, south_pole, andoyer);
  std::cout << "\tSemiminor distance:\t" << andoyer_minor << "\n";

  std::cout << "Vincenty function:\n";
  vincenty_major = boost::geometry::distance(equator_east, equator_west, vincenty);
  std::cout << "\tSemimajor distance:\t" << vincenty_major << "\n";
  vincenty_minor = boost::geometry::distance(north_pole, south_pole, vincenty);
  std::cout << "\tSemiminor distance:\t" << vincenty_minor << "\n";

  return 0;
}

Answer 1

作为替代方案，请查看 Charles FF Karney 的geoliblib。正如文档所述：“重点是返回准确的结果，误差接近四舍五入（约 5\xe2\x80\x9315 纳米）。”

\n