Modify ↓
Opened 17 months ago
Closed 3 months ago
#11931 closed Bugs (fixed)
spherical area wrong with pole encircling polygons, etc.
Reported by: | Charles Karney <charles@…> | Owned by: | barendgehrels |
---|---|---|---|
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:
- 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).
- 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 3 months ago by awulkiew
- Keywords area added
- Milestone changed from To Be Determined to Boost 1.64.0
- Resolution set to fixed
- Status changed from new to closed
Note: See
TracTickets for help on using
tickets.
Thanks!
Fixed by this PR: https://github.com/boostorg/geometry/pull/359
The new results are: