spherical area wrong with pole encircling polygons, etc.

[Charles Karney <charles@…>]

The following code illustrates a problem with handling of meridian crossing and pole encircling when computing spherical areas. It computes the area of a 4 identical diamond shaped areas at different positions on the sphere. Only the area of the last of these is computed correctly. The others are badly wrong.

#include <iostream>
#include <string>
#include <boost/geometry.hpp>

namespace bg = boost::geometry;
typedef bg::model::point<double, 2,
                         bg::cs::spherical_equatorial<bg::degree> > pt;

int main() {
  std::string polys[] = {
    " 0  89,  90  89,  180  89,  270  89",
    " 0 -89, -90 -89, -180 -89, -270 -89",
    "-1   0,   0  -1,    1   0,    0   1",
    "89   0,  90  -1,   91   0,   90   1"
  };
  // Correct area for all these polygons is
  double area0 = 0.00060926577032113655;
  bg::model::polygon<pt, false, false> poly;
  std::cout << "  area      relerr" << "\n";
  for (int i = 0; i < 4; ++i) {
    bg::read_wkt("POLYGON((" + polys[i] + "))", poly);
    double area = bg::area(poly);
    std::cout << area << " " << (area - area0) / area0 << "\n";
  }
}

Running this code gives

  area      relerr
0.000304633 -0.5
-6.28288 -10313.2
0 -1
0.000609266 2.31338e-15

Discussion:

• I can't make out the handling of this issue in the current code at all. It seems to be irredeemably broken.
• You need to do two things to handle these issues:
  1. Make sure that the sign of lambda_2 - lambda_1 correctly corresponds to east and west going edges (with some consistent convention of edges which pass over the pole).
  2. Count the number of times the polygon encircles a pole, e.g., by keeping track of when an edge crosses the prime meridian (again with some consistent treatment of cases when one or both ends of an edge lie on the prime meridian).
• If this count is odd, then the total area needs to be adjusted by 2*pi, half the area of the sphere.
• I believe I've handled these issues correctly in the PolygonAreaT class in ​GeographicLib.

[awulkiew]

Thanks!

Fixed by this PR: ​https://github.com/boostorg/geometry/pull/359

The new results are:

  area      relerr
0.000609266 -1.40992e-12
0.000609266 -1.40992e-12
0.000609266 1.77952e-16
0.000609266 1.42362e-15