Issues calculating geod.Inverse antipodal and near-antipodal distance

[Anonymous]

Upon raising an issue with the boost::geometry wgs84 distance functions on StackOverflow (see why is boost::geometry geographic vincenty distance inaccurate around the equator?
I was recommended to try geographiclib instead.

However, upon testing geographiclib with the same antipodal and near-antipodal points I got the following results:

GeographicLib near antipodal
    Semimajor distance:     20003931.458625
    Semiminor distance:     nan

GeographicLib antipodal
    Semimajor distance:     20003931.458625
    Semiminor distance:     0.000000

GeographicLib verify distance
    JFK to LHR distance:    5551759.400319

The Semimajor distance (i.e. around the equator) is inaccurate, by my calcualtions it should be 20037508.342789. And the Semiminor distance (i.e.from pole to pole) is a major issue.

Here's the code I used:

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

I'm sorry to be the bearer of bad news.

[Charles Karney]

You have the position of the south pole wrong. Please try applying this patch:

--- code.cpp.orig       2016-01-17 10:28:35.363515668 -0500
+++ code.cpp    2016-01-17 10:27:51.907070615 -0500
@@ -52,7 +52,7 @@
   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);
+  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);
@@ -67,7 +67,7 @@
   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);
+  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);
@@ -83,4 +83,3 @@

   return 0;
 }