Modify

Opened 22 months ago

Closed 8 months ago

#11931 closed Bugs (fixed)

spherical area wrong with pole encircling polygons, etc.

Reported by: Charles Karney <charles@…> Owned by: Barend Gehrels
Milestone: Boost 1.64.0 Component: geometry
Version: Boost 1.58.0 Severity: Showstopper
Keywords: area Cc: charles.karney@…

Description

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.

Attachments (0)

Change History (1)

comment:1 Changed 8 months ago by awulkiew

Keywords: area added
Milestone: To Be DeterminedBoost 1.64.0
Resolution: fixed
Status: newclosed

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

Modify Ticket

Change Properties
Set your email in Preferences
Action
as closed The owner will remain Barend Gehrels.
The resolution will be deleted.

Add Comment


E-mail address and name can be saved in the Preferences.

 
Note: See TracTickets for help on using tickets.