[Boost-bugs] [Boost C++ Libraries] #11931: spherical area wrong with pole encircling polygons, etc.

Subject: [Boost-bugs] [Boost C++ Libraries] #11931: spherical area wrong with pole encircling polygons, etc.
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2016-01-21 19:10:18

#11931: spherical area wrong with pole encircling polygons, etc.
----------------------------------------+---------------------------
 Reporter:  Charles Karney <charles@…>  |      Owner:  barendgehrels
     Type:  Bugs                        |     Status:  new
Milestone:  To Be Determined            |  Component:  geometry
  Version:  Boost 1.58.0                |   Severity:  Showstopper
 Keywords:                              |
----------------------------------------+---------------------------
 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.

 {{{#!c++
 #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 [http://geographiclib.sourceforge.net/html/ GeographicLib].

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/11931>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.