为什么 boost::geometry 赤道周围的文森蒂地理距离不准确？
Why is boost::geometry geographic Vincenty distance inaccurate around the Equator?
我需要一个函数来计算一对 WGS 84 positions to a high degree of accuracy and I was planning to use the geographic functions from boost geometry.
之间的距离
boost geometry Design Rational 状态：
There is the Andoyer method, fast and precise, and there is the Vincenty method, slower and more precise..
但是，当使用 Andoyer 和 Vincenty 策略测试 boost::geometry::distance 函数时，我得到了以下结果：
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 函数 returns 零！
问题出在：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;
}
作为替代方案，请查看 Charles F. F. Karney 的 geographiclib。正如文档所说："The emphasis is on returning accurate results with errors close to round-off (about 5–15 nanometers)."
我听从了@jwd630 的建议并查看了 geographiclib。
以下是结果：
WGS 84 values (metres)
Semimajor distance: 20037508.342789
Semiminor distance: 19970326.371123
GeographicLib near antipodal
Semimajor distance: 20003931.458625
Semiminor distance: 20003931.455275
GeographicLib antipodal
Semimajor distance: 20003931.458625
Semiminor distance: 20003931.458625
GeographicLib verify
JFK to LHR distance: 5551759.400319
即它为两极之间的半短距离（至 5dp）提供与 Vincenty 相同的距离，并计算赤道对映点的相同距离。
这是正确的，因为赤道对映点之间的最短距离是通过其中一个极点，而不是像默认提升 Andoyer 算法计算的那样围绕赤道。
所以@jwd630 上面的回答是正确的，在这三种算法中，geographiclib 是唯一一种计算整个 WGS84 大地水准面的正确距离的算法。
这里是测试代码：
/// GeographicLib WGS84 distance
// Note: M_PI is not part of the C or C++ standards, _USE_MATH_DEFINES enables it
#define _USE_MATH_DEFINES
#include <GeographicLib/Geodesic.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*/)
{
const GeographicLib::Geodesic& geod(GeographicLib::Geodesic::WGS84());
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 boost Andoyer to fail.
const double DELTA(0.000000015);
std::pair<double, double> near_equator_east (0.0, 90.0 - DELTA);
std::pair<double, double> near_equator_west (0.0, -90.0 + DELTA);
std::cout << "GeographicLib near antipodal\n";
double distance_metres(0.0);
geod.Inverse(near_equator_east.first, near_equator_east.second,
near_equator_west.first, near_equator_west.second, distance_metres);
std::cout << "\tSemimajor distance:\t" << distance_metres << "\n";
std::pair<double, double> near_north_pole (90.0 - DELTA, 0.0);
std::pair<double, double> near_south_pole (-90.0 + DELTA, 0.0);
geod.Inverse(near_north_pole.first, near_north_pole.second,
near_south_pole.first, near_south_pole.second, distance_metres);
std::cout << "\tSemiminor distance:\t" << distance_metres << "\n\n";
std::pair<double, double> equator_east (0.0, 90.0);
std::pair<double, double> equator_west (0.0, -90.0);
std::cout << "GeographicLib antipodal\n";
geod.Inverse(equator_east.first, equator_east.second,
equator_west.first, equator_west.second, distance_metres);
std::cout << "\tSemimajor distance:\t" << distance_metres << "\n";
std::pair<double, double> north_pole (90.0, 0.0);
std::pair<double, double> south_pole (-90.0, 0.0);
geod.Inverse(north_pole.first, north_pole.second,
south_pole.first, south_pole.second, distance_metres);
std::cout << "\tSemiminor distance:\t" << distance_metres << "\n\n";
std::pair<double, double> JFK (40.6, -73.8);
std::pair<double, double> LHR (51.6, -0.5);
std::cout << "GeographicLib verify distance\n";
geod.Inverse(JFK.first, JFK.second,
LHR.first, LHR.second, distance_metres);
std::cout << "\tJFK to LHR distance:\t" << distance_metres << std::endl;
return 0;
}
在他的论文中Algorithms for geodesics,
Charles F. F. Karney 指出 "Vincenty’s method fails to converge for nearly antipodal points"。
这可能会回答我原来的问题，即 Vincenty 算法不适合对映点。
注意：我提出了 boost 票 #11817 描述了问题所在
Andoyer 算法 returns 对映点为零，并向 boost 发送了一个拉取请求并对其进行了修复。
但是，对于错误距离的唯一正确修复是使用正确的算法，即：geographiclib
非常感谢 Charles F. F. Karney (@cffk) 礼貌地指出我的愚蠢错误！
我需要一个函数来计算一对 WGS 84 positions to a high degree of accuracy and I was planning to use the geographic functions from boost geometry.
boost geometry Design Rational 状态：
There is the Andoyer method, fast and precise, and there is the Vincenty method, slower and more precise..
但是，当使用 Andoyer 和 Vincenty 策略测试 boost::geometry::distance 函数时，我得到了以下结果：
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 函数 returns 零！
问题出在：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;
}
作为替代方案，请查看 Charles F. F. Karney 的 geographiclib。正如文档所说："The emphasis is on returning accurate results with errors close to round-off (about 5–15 nanometers)."
我听从了@jwd630 的建议并查看了 geographiclib。
以下是结果：
WGS 84 values (metres)
Semimajor distance: 20037508.342789
Semiminor distance: 19970326.371123
GeographicLib near antipodal
Semimajor distance: 20003931.458625
Semiminor distance: 20003931.455275
GeographicLib antipodal
Semimajor distance: 20003931.458625
Semiminor distance: 20003931.458625
GeographicLib verify
JFK to LHR distance: 5551759.400319
即它为两极之间的半短距离（至 5dp）提供与 Vincenty 相同的距离，并计算赤道对映点的相同距离。
这是正确的，因为赤道对映点之间的最短距离是通过其中一个极点，而不是像默认提升 Andoyer 算法计算的那样围绕赤道。
所以@jwd630 上面的回答是正确的，在这三种算法中，geographiclib 是唯一一种计算整个 WGS84 大地水准面的正确距离的算法。
这里是测试代码：
/// GeographicLib WGS84 distance
// Note: M_PI is not part of the C or C++ standards, _USE_MATH_DEFINES enables it
#define _USE_MATH_DEFINES
#include <GeographicLib/Geodesic.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*/)
{
const GeographicLib::Geodesic& geod(GeographicLib::Geodesic::WGS84());
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 boost Andoyer to fail.
const double DELTA(0.000000015);
std::pair<double, double> near_equator_east (0.0, 90.0 - DELTA);
std::pair<double, double> near_equator_west (0.0, -90.0 + DELTA);
std::cout << "GeographicLib near antipodal\n";
double distance_metres(0.0);
geod.Inverse(near_equator_east.first, near_equator_east.second,
near_equator_west.first, near_equator_west.second, distance_metres);
std::cout << "\tSemimajor distance:\t" << distance_metres << "\n";
std::pair<double, double> near_north_pole (90.0 - DELTA, 0.0);
std::pair<double, double> near_south_pole (-90.0 + DELTA, 0.0);
geod.Inverse(near_north_pole.first, near_north_pole.second,
near_south_pole.first, near_south_pole.second, distance_metres);
std::cout << "\tSemiminor distance:\t" << distance_metres << "\n\n";
std::pair<double, double> equator_east (0.0, 90.0);
std::pair<double, double> equator_west (0.0, -90.0);
std::cout << "GeographicLib antipodal\n";
geod.Inverse(equator_east.first, equator_east.second,
equator_west.first, equator_west.second, distance_metres);
std::cout << "\tSemimajor distance:\t" << distance_metres << "\n";
std::pair<double, double> north_pole (90.0, 0.0);
std::pair<double, double> south_pole (-90.0, 0.0);
geod.Inverse(north_pole.first, north_pole.second,
south_pole.first, south_pole.second, distance_metres);
std::cout << "\tSemiminor distance:\t" << distance_metres << "\n\n";
std::pair<double, double> JFK (40.6, -73.8);
std::pair<double, double> LHR (51.6, -0.5);
std::cout << "GeographicLib verify distance\n";
geod.Inverse(JFK.first, JFK.second,
LHR.first, LHR.second, distance_metres);
std::cout << "\tJFK to LHR distance:\t" << distance_metres << std::endl;
return 0;
}
在他的论文中Algorithms for geodesics,
Charles F. F. Karney 指出 "Vincenty’s method fails to converge for nearly antipodal points"。
这可能会回答我原来的问题，即 Vincenty 算法不适合对映点。
注意：我提出了 boost 票 #11817 描述了问题所在
Andoyer 算法 returns 对映点为零，并向 boost 发送了一个拉取请求并对其进行了修复。
但是，对于错误距离的唯一正确修复是使用正确的算法，即：geographiclib
非常感谢 Charles F. F. Karney (@cffk) 礼貌地指出我的愚蠢错误！