GeographicLib  2.1.2
PolygonArea.hpp
Go to the documentation of this file.
1 /**
2  * \file PolygonArea.hpp
3  * \brief Header for GeographicLib::PolygonAreaT class
4  *
5  * Copyright (c) Charles Karney (2010-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_POLYGONAREA_HPP)
11 #define GEOGRAPHICLIB_POLYGONAREA_HPP 1
12 
15 #include <GeographicLib/Rhumb.hpp>
17 
18 namespace GeographicLib {
19 
20  /**
21  * \brief Polygon areas
22  *
23  * This computes the area of a polygon whose edges are geodesics using the
24  * method given in Section 6 of
25  * - C. F. F. Karney,
26  * <a href="https://doi.org/10.1007/s00190-012-0578-z">
27  * Algorithms for geodesics</a>,
28  * J. Geodesy <b>87</b>, 43--55 (2013);
29  * DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
30  * 10.1007/s00190-012-0578-z</a>;
31  * addenda:
32  * <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
33  * geod-addenda.html</a>.
34  *
35  * Arbitrarily complex polygons are allowed. In the case self-intersecting
36  * of polygons the area is accumulated "algebraically", e.g., the areas of
37  * the 2 loops in a figure-8 polygon will partially cancel.
38  *
39  * This class lets you add vertices and edges one at a time to the polygon.
40  * The sequence must start with a vertex and thereafter vertices and edges
41  * can be added in any order. Any vertex after the first creates a new edge
42  * which is the \e shortest geodesic from the previous vertex. In some
43  * cases there may be two or many such shortest geodesics and the area is
44  * then not uniquely defined. In this case, either add an intermediate
45  * vertex or add the edge \e as an edge (by defining its direction and
46  * length).
47  *
48  * The area and perimeter are accumulated at two times the standard floating
49  * point precision to guard against the loss of accuracy with many-sided
50  * polygons. At any point you can ask for the perimeter and area so far.
51  * There's an option to treat the points as defining a polyline instead of a
52  * polygon; in that case, only the perimeter is computed.
53  *
54  * This is a templated class to allow it to be used with Geodesic,
55  * GeodesicExact, and Rhumb. GeographicLib::PolygonArea,
56  * GeographicLib::PolygonAreaExact, and GeographicLib::PolygonAreaRhumb are
57  * typedefs for these cases.
58  *
59  * For GeographicLib::PolygonArea (edges defined by Geodesic), an upper bound
60  * on the error is about 0.1 m<sup>2</sup> per vertex. However this is a
61  * wildly pessimistic estimate in most cases. A more realistic estimate of
62  * the error is given by a test involving 10<sup>7</sup> approximately
63  * regular polygons on the WGS84 ellipsoid. The centers and the orientations
64  * of the polygons were uniformly distributed, the number of vertices was
65  * log-uniformly distributed in [3, 300], and the center to vertex distance
66  * log-uniformly distributed in [0.1 m, 9000 km].
67  *
68  * Using double precision (the standard precision for GeographicLib), the
69  * maximum error in the perimeter was 200 nm, and the maximum error in the
70  * area was<pre>
71  * 0.0013 m^2 for perimeter < 10 km
72  * 0.0070 m^2 for perimeter < 100 km
73  * 0.070 m^2 for perimeter < 1000 km
74  * 0.11 m^2 for all perimeters
75  * </pre>
76  * The errors are given in terms of the perimeter, because it is expected
77  * that the errors depend mainly on the number of edges and the edge lengths.
78  *
79  * Using long doubles (GEOGRPAHICLIB_PRECISION = 3), the maximum error in the
80  * perimeter was 200 pm, and the maximum error in the area was<pre>
81  * 0.7 mm^2 for perim < 10 km
82  * 3.2 mm^2 for perimeter < 100 km
83  * 21 mm^2 for perimeter < 1000 km
84  * 45 mm^2 for all perimeters
85  * </pre>
86  *
87  * @tparam GeodType the geodesic class to use.
88  *
89  * Example of use:
90  * \include example-PolygonArea.cpp
91  *
92  * <a href="Planimeter.1.html">Planimeter</a> is a command-line utility
93  * providing access to the functionality of PolygonAreaT.
94  **********************************************************************/
95 
96  template<class GeodType = Geodesic>
97  class PolygonAreaT {
98  private:
99  typedef Math::real real;
100  GeodType _earth;
101  real _area0; // Full ellipsoid area
102  bool _polyline; // Assume polyline (don't close and skip area)
103  unsigned _mask;
104  unsigned _num;
105  int _crossings;
106  Accumulator<> _areasum, _perimetersum;
107  real _lat0, _lon0, _lat1, _lon1;
108  static int transit(real lon1, real lon2);
109  // an alternate version of transit to deal with longitudes in the direct
110  // problem.
111  static int transitdirect(real lon1, real lon2);
112  void Remainder(Accumulator<>& a) const { a.remainder(_area0); }
113  void Remainder(real& a) const {
114  using std::remainder;
115  a = remainder(a, _area0);
116  }
117  template<typename T>
118  void AreaReduce(T& area, int crossings, bool reverse, bool sign) const;
119  public:
120 
121  /**
122  * Constructor for PolygonAreaT.
123  *
124  * @param[in] earth the Geodesic object to use for geodesic calculations.
125  * @param[in] polyline if true that treat the points as defining a polyline
126  * instead of a polygon (default = false).
127  **********************************************************************/
128  PolygonAreaT(const GeodType& earth, bool polyline = false)
129  : _earth(earth)
130  , _area0(_earth.EllipsoidArea())
131  , _polyline(polyline)
132  , _mask(GeodType::LATITUDE | GeodType::LONGITUDE | GeodType::DISTANCE |
133  (_polyline ? GeodType::NONE :
134  GeodType::AREA | GeodType::LONG_UNROLL))
135  { Clear(); }
136 
137  /**
138  * Clear PolygonAreaT, allowing a new polygon to be started.
139  **********************************************************************/
140  void Clear() {
141  _num = 0;
142  _crossings = 0;
143  _areasum = 0;
144  _perimetersum = 0;
145  _lat0 = _lon0 = _lat1 = _lon1 = Math::NaN();
146  }
147 
148  /**
149  * Add a point to the polygon or polyline.
150  *
151  * @param[in] lat the latitude of the point (degrees).
152  * @param[in] lon the longitude of the point (degrees).
153  *
154  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
155  **********************************************************************/
156  void AddPoint(real lat, real lon);
157 
158  /**
159  * Add an edge to the polygon or polyline.
160  *
161  * @param[in] azi azimuth at current point (degrees).
162  * @param[in] s distance from current point to next point (meters).
163  *
164  * This does nothing if no points have been added yet. Use
165  * PolygonAreaT::CurrentPoint to determine the position of the new vertex.
166  **********************************************************************/
167  void AddEdge(real azi, real s);
168 
169  /**
170  * Return the results so far.
171  *
172  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
173  * traversal counts as a positive area.
174  * @param[in] sign if true then return a signed result for the area if
175  * the polygon is traversed in the "wrong" direction instead of returning
176  * the area for the rest of the earth.
177  * @param[out] perimeter the perimeter of the polygon or length of the
178  * polyline (meters).
179  * @param[out] area the area of the polygon (meters<sup>2</sup>); only set
180  * if \e polyline is false in the constructor.
181  * @return the number of points.
182  *
183  * More points can be added to the polygon after this call.
184  **********************************************************************/
185  unsigned Compute(bool reverse, bool sign,
186  real& perimeter, real& area) const;
187 
188  /**
189  * Return the results assuming a tentative final test point is added;
190  * however, the data for the test point is not saved. This lets you report
191  * a running result for the perimeter and area as the user moves the mouse
192  * cursor. Ordinary floating point arithmetic is used to accumulate the
193  * data for the test point; thus the area and perimeter returned are less
194  * accurate than if PolygonAreaT::AddPoint and PolygonAreaT::Compute are
195  * used.
196  *
197  * @param[in] lat the latitude of the test point (degrees).
198  * @param[in] lon the longitude of the test point (degrees).
199  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
200  * traversal counts as a positive area.
201  * @param[in] sign if true then return a signed result for the area if
202  * the polygon is traversed in the "wrong" direction instead of returning
203  * the area for the rest of the earth.
204  * @param[out] perimeter the approximate perimeter of the polygon or length
205  * of the polyline (meters).
206  * @param[out] area the approximate area of the polygon
207  * (meters<sup>2</sup>); only set if polyline is false in the
208  * constructor.
209  * @return the number of points.
210  *
211  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
212  **********************************************************************/
213  unsigned TestPoint(real lat, real lon, bool reverse, bool sign,
214  real& perimeter, real& area) const;
215 
216  /**
217  * Return the results assuming a tentative final test point is added via an
218  * azimuth and distance; however, the data for the test point is not saved.
219  * This lets you report a running result for the perimeter and area as the
220  * user moves the mouse cursor. Ordinary floating point arithmetic is used
221  * to accumulate the data for the test point; thus the area and perimeter
222  * returned are less accurate than if PolygonAreaT::AddEdge and
223  * PolygonAreaT::Compute are used.
224  *
225  * @param[in] azi azimuth at current point (degrees).
226  * @param[in] s distance from current point to final test point (meters).
227  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
228  * traversal counts as a positive area.
229  * @param[in] sign if true then return a signed result for the area if
230  * the polygon is traversed in the "wrong" direction instead of returning
231  * the area for the rest of the earth.
232  * @param[out] perimeter the approximate perimeter of the polygon or length
233  * of the polyline (meters).
234  * @param[out] area the approximate area of the polygon
235  * (meters<sup>2</sup>); only set if polyline is false in the
236  * constructor.
237  * @return the number of points.
238  **********************************************************************/
239  unsigned TestEdge(real azi, real s, bool reverse, bool sign,
240  real& perimeter, real& area) const;
241 
242  /** \name Inspector functions
243  **********************************************************************/
244  ///@{
245  /**
246  * @return \e a the equatorial radius of the ellipsoid (meters). This is
247  * the value inherited from the Geodesic object used in the constructor.
248  **********************************************************************/
249 
250  Math::real EquatorialRadius() const { return _earth.EquatorialRadius(); }
251 
252  /**
253  * @return \e f the flattening of the ellipsoid. This is the value
254  * inherited from the Geodesic object used in the constructor.
255  **********************************************************************/
256  Math::real Flattening() const { return _earth.Flattening(); }
257 
258  /**
259  * Report the previous vertex added to the polygon or polyline.
260  *
261  * @param[out] lat the latitude of the point (degrees).
262  * @param[out] lon the longitude of the point (degrees).
263  *
264  * If no points have been added, then NaNs are returned. Otherwise, \e lon
265  * will be in the range [&minus;180&deg;, 180&deg;].
266  **********************************************************************/
267  void CurrentPoint(real& lat, real& lon) const
268  { lat = _lat1; lon = _lon1; }
269 
270  /**
271  * Report the number of points currently in the polygon or polyline.
272  *
273  * @return the number of points.
274  *
275  * If no points have been added, then 0 is returned.
276  **********************************************************************/
277  unsigned NumberPoints() const { return _num; }
278 
279  /**
280  * Report whether the current object is a polygon or a polyline.
281  *
282  * @return true if the object is a polyline.
283  **********************************************************************/
284  bool Polyline() const { return _polyline; }
285  ///@}
286  };
287 
288  /**
289  * @relates PolygonAreaT
290  *
291  * Polygon areas using Geodesic. This should be used if the flattening is
292  * small.
293  **********************************************************************/
295 
296  /**
297  * @relates PolygonAreaT
298  *
299  * Polygon areas using GeodesicExact. (But note that the implementation of
300  * areas in GeodesicExact uses a high order series and this is only accurate
301  * for modest flattenings.)
302  **********************************************************************/
304 
305  /**
306  * @relates PolygonAreaT
307  *
308  * Polygon areas using Rhumb.
309  **********************************************************************/
311 
312 } // namespace GeographicLib
313 
314 #endif // GEOGRAPHICLIB_POLYGONAREA_HPP
Header for GeographicLib::Accumulator class.
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Header for GeographicLib::GeodesicExact class.
Header for GeographicLib::Geodesic class.
Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes.
An accumulator for sums.
Definition: Accumulator.hpp:40
Accumulator & remainder(T y)
static T NaN()
Definition: Math.cpp:250
PolygonAreaT(const GeodType &earth, bool polyline=false)
void AddPoint(real lat, real lon)
Definition: PolygonArea.cpp:70
unsigned Compute(bool reverse, bool sign, real &perimeter, real &area) const
PolygonAreaT< Rhumb > PolygonAreaRhumb
void CurrentPoint(real &lat, real &lon) const
Math::real Flattening() const
Math::real EquatorialRadius() const
PolygonAreaT< Geodesic > PolygonArea
unsigned TestPoint(real lat, real lon, bool reverse, bool sign, real &perimeter, real &area) const
unsigned NumberPoints() const
void AddEdge(real azi, real s)
Definition: PolygonArea.cpp:89
PolygonAreaT< GeodesicExact > PolygonAreaExact
unsigned TestEdge(real azi, real s, bool reverse, bool sign, real &perimeter, real &area) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12