Why is boost::geometry geographic Vincenty distance inaccurate around the Equator?

Question

I require a function to calculate the distance between a pair of WGS 84 positions to a high degree of accuracy and I was planning to use the geographic functions from boost geometry.

The boost geometry Design Rational states:

There is the Andoyer method, fast and precise, and there is the Vincenty method, slower and more precise..

However, when testing the boost::geometry::distance function with both the Andoyer and Vincenty strategies, I got the following results:

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

The Vincenty distances along the Semimajor axis (i.e. around the Equator) are less than the distance around the Semiminor axis between the North and South Poles. That can't be correct.

The Semiminor and Andoyer distances look reasonable. Except when the points are on opposing side of the Earth, when the boost Andoyer function returns zero!

Is the problem in: the Vincenty algorithm, the boost geometry implementation of it, or my test code?

Test code:

/// 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

I followed the advice of @jwd630 and checked out geographiclib.
Here are the results:

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

I.e. it provides the same distance as Vincenty for the Semiminor distance between the poles (to 5dp) and it calculates the same distance for antipodal points at the Equator.

This is correct, since the shortest distance between the antipodal points at the Equator is via one of the Poles, not around the equator as the default boost Andoyer algorithm calculates.

So @jwd630's answer above is right and of the three algorithms, geographiclib is the only one to calculate the correct distance over the whole WGS84 geoid.

Here is the test code:

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

In his paper Algorithms for geodesics, Charles F. F. Karney states that "Vincenty’s method fails to converge for nearly antipodal points". Which may answer my original question, i.e. that the Vincenty algorithm isn't suitable for antipodal points.

Note: I've raised boost ticket #11817 describing the issue where the Andoyer algorithm returns zero for antipodal points and sent a pull request to boost with a fix for it.

However, the only correct fix for the incorrect distances is to use the correct algorithm, namely: geographiclib

Many thanks to Charles F. F. Karney (@cffk) for politely pointing out my silly mistakes!

Answer 2

As an alternative check out Charles F. F. Karney's geographiclib. As the documentation says: "The emphasis is on returning accurate results with errors close to round-off (about 5–15 nanometers)."