Почему boost :: геометрия географического Vincenty расстояние неточно вокруг экватора?

Question

Мне нужна функция для вычисления расстояния между двумя позициями WGS 84 с высокой степенью точности, и я планировал использовать функции geographic от boost geometry.

В boost geometry Design Rational состояния:

Там является метод 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 вдоль оси полуосновника (т. вокруг экватора) меньше расстояния вокруг оси Семиминор между северным и южным полюсами. Это не может быть правильно.

Полуминор и Andoyer расстояния выглядят разумно. За исключением случаев, когда точки находятся на противоположной стороне Земли, когда функция boostAndoyer возвращает ноль!

Является ли проблема в: алгоритме 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

В качестве альтернативного заканчивало Чарльз Ф. Ф. Karney-х geographiclib. Как говорится в документации: «Акцент делается на возвращении точных результатов с ошибками, близкими к округлению (около 5-15 нанометров)».

Answer 2

Я последовал совету @ 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), и он вычисляет такое же расстояние для антиподальных точек на экваторе.

Это правильно, так как кратчайшее расстояние между антиподальными точками на экваторе осуществляется через один из полюсов, а не вокруг экватора, поскольку алгоритм вычисления по умолчанию 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, Чарльз Ф. Ф. Карни утверждает, что «метод Vincenty в не сходится почти диаметрально противоположных точек». Что может ответить на мой первоначальный вопрос, т. Е. Что алгоритм Vincenty не подходит для антиподальных точек.

Примечания: Я поднял boost билет #11817 описывающий вопрос, где в Andoyer алгоритм возвращает ноль для антипода точек и послал запрос тянуть к boost с исправлением для него.

Однако единственно правильное исправление неправильных расстояний использовать правильный алгоритм, а именно: geographiclib

Большое спасибо Чарльзу Ф. Ф. Карни (@cffk) для вежливо указывая на мои глупые ошибки!