karney_direct longitude is inaccurate

[kenba]

The karney_direct formula (formulas/karney_direct.hpp) is an implementation of Charles Karney’s GeographicLib Geodesic::Direct function. Accuracy is a key feature of the GeographicLib functions, see #449.

I therefore expect the karney_direct formula to produce results: azimuth, latitude and longitude, that are very close to those produced by the corresponding GeographicLib Geodesic::Direct function.

Comparing the output of the karney_direct formula with Geodesic::Direct over a grid of latitudes and longitudes, the reverse_azimuths and latitudes match the expected values quite closely but the longitudes are less accurate by over 4 orders of magnitude!, see program below.

I believe that the longitude inaccuracy may be caused by a bug in the evaluate_coeffs_C3 function in the util/series_expansion.hpp file: the variable m is not being passed to the math::polyval function.

The comment above "evaluate_coeffs_C3x" states: "The C++ code ... is generated by the following Maxima script..." Therefore, this issue may be caused by Maxima or the script.

#include <boost/geometry/algorithms/detail/azimuth.hpp>
#include <boost/geometry/formulas/karney_direct.hpp>
#include <GeographicLib/Geodesic.hpp>
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <iostream>

namespace
{
  constexpr long double wgs84_a{ 6378137.0L };
  constexpr long double wgs84_f{ 1.0L / 298.257223563L };
  constexpr long double one_minus_f{ 1.0L - wgs84_f };
  constexpr long double wgs84_b{ wgs84_a * one_minus_f };

  const boost::geometry::srs::spheroid<long double> BoostWGS84(wgs84_a, wgs84_b);

  // boost karney_direct function class with azimuth output and SeriesOrder = 6
  template <typename T>
  using BoostKarneyDirect_6 = boost::geometry::formula::karney_direct <T, true, true, false, false, 6u>;

  // boost karney_direct function class with azimuth output and SeriesOrder = 8
  template <typename T>
  using BoostKarneyDirect_8 = boost::geometry::formula::karney_direct <T, true, true, false, false, 8u>;

  // boost test BOOST_CHECK_CLOSE macro takes a percentage accuracy parameter
  const auto EPSILON(std::numeric_limits<double>::epsilon());
  const auto CALCULATION_TOLERANCE(100 * EPSILON);
}

BOOST_AUTO_TEST_SUITE(Test_boost_karney_direct)

BOOST_AUTO_TEST_CASE(test_Geodesic_positions)
{
  const auto GeographicLibWGS84(GeographicLib::Geodesic::WGS84());

  // Loop around latitudes: 0 to 89 degrees
  for (auto i(0); i < 90; ++i)
  {
    // The latitude in degrees.
    double latitude(1.0 * i);

    // Loop around longitudes: 1 to 179 degrees
    for (auto j(1); j < 180; ++j)
    {
      // The longitude in degrees.
      double longitude(1.0 * j);

      // The Geodesic: distance in metres, start azimuth and finish azimuth in degrees.
      double distance_m, azimuth, azi2;
      GeographicLibWGS84.Inverse(0.0, 0.0, latitude, longitude, distance_m, azimuth, azi2);

      // The GeographicLib position and azimuth at the distance in metres
      double lat2k, lon2k, azi2k;
      GeographicLibWGS84.Direct(0.0, 0.0, azimuth, distance_m, lat2k, lon2k, azi2k);
      BOOST_CHECK_CLOSE(latitude, lat2k, 140 * CALCULATION_TOLERANCE);
      BOOST_CHECK_CLOSE(longitude, lon2k, 120 * CALCULATION_TOLERANCE);

      // The boost karney_direct order 6 position at the azimuth and distance in metres.
      auto results_6(BoostKarneyDirect_6<double>::apply(0.0, 0.0, distance_m, azimuth, BoostWGS84));
      BOOST_CHECK_CLOSE(azi2, results_6.reverse_azimuth, 140 * CALCULATION_TOLERANCE);
      BOOST_CHECK_CLOSE(latitude, results_6.lat2, 220 * CALCULATION_TOLERANCE);

      /******** Test below only passes with >= 10172000 * CALCULATION_TOLERANCE !! ********/
      BOOST_CHECK_CLOSE(longitude, results_6.lon2, 10171000 * CALCULATION_TOLERANCE);
      /*****************************************************************************/

      // The boost karney_direct order 8 position at the azimuth and distance in metres.
      auto results_8(BoostKarneyDirect_8<double>::apply(0.0, 0.0, distance_m, azimuth, BoostWGS84));
      BOOST_CHECK_CLOSE(azi2, results_8.reverse_azimuth, 140 * CALCULATION_TOLERANCE);
      BOOST_CHECK_CLOSE(latitude, results_8.lat2, 220 * CALCULATION_TOLERANCE);

      /******** Test below only passes with >= 10174000 * CALCULATION_TOLERANCE !! ********/
      BOOST_CHECK_CLOSE(longitude, results_8.lon2, 10173000 * CALCULATION_TOLERANCE);
      /*****************************************************************************/
    }
  }
}

BOOST_AUTO_TEST_SUITE_END()

[cffk]

Indeed it looks like the call to math::polyval on line 645 of geometry/util/series_expansion.hpp is wrong. I think

math::polyval(coeffs2.begin(), coeffs2.begin() + offset, eps);

should be replaced by

math::polyval(coeffs2.begin() + offset,
              coeffs2.begin() + offset + m + 1,
              eps);

and maybe offset needs to be initialized to 0 instead of 1 in this function.

[kenba]

I attempted to implement the changes proposed by @cffk above yesterday (including initialising offset to 0), only to discover that there is also an issue with polyval and the polynomial coefficients.

math::polyval evaluates polynomial coefficients in descending eps order, while evaluate_coeffs_C3x (and the other functions in series_expansion.hpp) define their polynomial coefficients in ascending eps order.

Having the polynomial coefficients in ascending eps order makes them easier to compare with the values in @cffk's paper: Geodesics on an ellipsoid of revolution. However, the GeographicLib Geodesic.cpp file has the coefficients descending eps order, so it is easier to compare these values in descending eps order.

I compared some of the polynomial coefficients in evaluate_coeffs_C3 with those in @cffk's paper and Geodesic.cpp and found that evaluate_coeffs_C3 includes extra coefficients...

This issue could be resolved either by changing math::polyval and/or reordering/recreating the polynomial coefficients in addition to the above changes.

Note: math::horner_evaluate evaluates polynomial coefficients in ascending eps order, so if the polynomial coefficients were changed to descending order, then this function would also need to be changed.

[kenba]

I have just committed a pull request to fix this issue.

It contains the changes proposed by @cffk and calls horner_evaluate instead of polyval to evaluate the polynomial coefficients in ascending eps order, i.e. it is the minimum change to fix the issue.

With this change, the example tests pass with both longitude tolerances set to 130 * CALCULATION_TOLERANCE.

[vissarion]

https://github.com/boostorg/geometry/pull/562 is merged. Should I close the issue?

[kenba]

I believe so, but then again I wrote the fix! Charles @cffk is the expert, I would ask his permission if I were you.

[cffk]

I haven't been paying much attention to the boost geometry package lately. So I'm not in a position to provide a decent response here. However, I can recommend a piece of due diligence...

The fact that the original boost implementation goofed on evaluating the series expansion for longitude suggests that the other series expansions used in the solution of the goedesic problem be scrutinized for similar problems too. The other expansions are those for the distance (and the corresponding reverted series), the reduced length + geodesic scale, and the series for the area. (I'm not sure if all of these are used by boost geometry.)

[vissarion]

It is not a priority but indeed this (accuracy of geodesic algorithms in Boost Geometry) is something that should be tested more systematically. Thanks for the feedback. I am closing the issue.