GeographicLib  1.50.1
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-2019) <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  * @tparam GeodType the geodesic class to use.
60  *
61  * Example of use:
62  * \include example-PolygonArea.cpp
63  *
64  * <a href="Planimeter.1.html">Planimeter</a> is a command-line utility
65  * providing access to the functionality of PolygonAreaT.
66  **********************************************************************/
67 
68  template <class GeodType = Geodesic>
69  class PolygonAreaT {
70  private:
71  typedef Math::real real;
72  GeodType _earth;
73  real _area0; // Full ellipsoid area
74  bool _polyline; // Assume polyline (don't close and skip area)
75  unsigned _mask;
76  unsigned _num;
77  int _crossings;
78  Accumulator<> _areasum, _perimetersum;
79  real _lat0, _lon0, _lat1, _lon1;
80  static int transit(real lon1, real lon2) {
81  // Return 1 or -1 if crossing prime meridian in east or west direction.
82  // Otherwise return zero.
83  // Compute lon12 the same way as Geodesic::Inverse.
84  lon1 = Math::AngNormalize(lon1);
85  lon2 = Math::AngNormalize(lon2);
86  real lon12 = Math::AngDiff(lon1, lon2);
87  // Treat 0 as negative in these tests. This balances +/- 180 being
88  // treated as positive, i.e., +180.
89  int cross =
90  lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
91  (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
92  return cross;
93  }
94  // an alternate version of transit to deal with longitudes in the direct
95  // problem.
96  static int transitdirect(real lon1, real lon2) {
97  // Compute exactly the parity of
98  // int(ceil(lon2 / 360)) - int(ceil(lon1 / 360))
99  lon1 = Math::remainder(lon1, real(720));
100  lon2 = Math::remainder(lon2, real(720));
101  return ( (lon2 <= 0 && lon2 > -360 ? 1 : 0) -
102  (lon1 <= 0 && lon1 > -360 ? 1 : 0) );
103  }
104  void Remainder(Accumulator<>& a) const { a.remainder(_area0); }
105  void Remainder(real& a) const { a = Math::remainder(a, _area0); }
106  template <typename T>
107  void AreaReduce(T& area, int crossings, bool reverse, bool sign) const;
108  public:
109 
110  /**
111  * Constructor for PolygonAreaT.
112  *
113  * @param[in] earth the Geodesic object to use for geodesic calculations.
114  * @param[in] polyline if true that treat the points as defining a polyline
115  * instead of a polygon (default = false).
116  **********************************************************************/
117  PolygonAreaT(const GeodType& earth, bool polyline = false)
118  : _earth(earth)
119  , _area0(_earth.EllipsoidArea())
120  , _polyline(polyline)
121  , _mask(GeodType::LATITUDE | GeodType::LONGITUDE | GeodType::DISTANCE |
122  (_polyline ? GeodType::NONE :
123  GeodType::AREA | GeodType::LONG_UNROLL))
124  { Clear(); }
125 
126  /**
127  * Clear PolygonAreaT, allowing a new polygon to be started.
128  **********************************************************************/
129  void Clear() {
130  _num = 0;
131  _crossings = 0;
132  _areasum = 0;
133  _perimetersum = 0;
134  _lat0 = _lon0 = _lat1 = _lon1 = Math::NaN();
135  }
136 
137  /**
138  * Add a point to the polygon or polyline.
139  *
140  * @param[in] lat the latitude of the point (degrees).
141  * @param[in] lon the longitude of the point (degrees).
142  *
143  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
144  **********************************************************************/
145  void AddPoint(real lat, real lon);
146 
147  /**
148  * Add an edge to the polygon or polyline.
149  *
150  * @param[in] azi azimuth at current point (degrees).
151  * @param[in] s distance from current point to next point (meters).
152  *
153  * This does nothing if no points have been added yet. Use
154  * PolygonAreaT::CurrentPoint to determine the position of the new vertex.
155  **********************************************************************/
156  void AddEdge(real azi, real s);
157 
158  /**
159  * Return the results so far.
160  *
161  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
162  * traversal counts as a positive area.
163  * @param[in] sign if true then return a signed result for the area if
164  * the polygon is traversed in the "wrong" direction instead of returning
165  * the area for the rest of the earth.
166  * @param[out] perimeter the perimeter of the polygon or length of the
167  * polyline (meters).
168  * @param[out] area the area of the polygon (meters<sup>2</sup>); only set
169  * if \e polyline is false in the constructor.
170  * @return the number of points.
171  *
172  * More points can be added to the polygon after this call.
173  **********************************************************************/
174  unsigned Compute(bool reverse, bool sign,
175  real& perimeter, real& area) const;
176 
177  /**
178  * Return the results assuming a tentative final test point is added;
179  * however, the data for the test point is not saved. This lets you report
180  * a running result for the perimeter and area as the user moves the mouse
181  * cursor. Ordinary floating point arithmetic is used to accumulate the
182  * data for the test point; thus the area and perimeter returned are less
183  * accurate than if PolygonAreaT::AddPoint and PolygonAreaT::Compute are
184  * used.
185  *
186  * @param[in] lat the latitude of the test point (degrees).
187  * @param[in] lon the longitude of the test point (degrees).
188  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
189  * traversal counts as a positive area.
190  * @param[in] sign if true then return a signed result for the area if
191  * the polygon is traversed in the "wrong" direction instead of returning
192  * the area for the rest of the earth.
193  * @param[out] perimeter the approximate perimeter of the polygon or length
194  * of the polyline (meters).
195  * @param[out] area the approximate area of the polygon
196  * (meters<sup>2</sup>); only set if polyline is false in the
197  * constructor.
198  * @return the number of points.
199  *
200  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
201  **********************************************************************/
202  unsigned TestPoint(real lat, real lon, bool reverse, bool sign,
203  real& perimeter, real& area) const;
204 
205  /**
206  * Return the results assuming a tentative final test point is added via an
207  * azimuth and distance; however, the data for the test point is not saved.
208  * This lets you report a running result for the perimeter and area as the
209  * user moves the mouse cursor. Ordinary floating point arithmetic is used
210  * to accumulate the data for the test point; thus the area and perimeter
211  * returned are less accurate than if PolygonAreaT::AddEdge and
212  * PolygonAreaT::Compute are used.
213  *
214  * @param[in] azi azimuth at current point (degrees).
215  * @param[in] s distance from current point to final test point (meters).
216  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
217  * traversal counts as a positive area.
218  * @param[in] sign if true then return a signed result for the area if
219  * the polygon is traversed in the "wrong" direction instead of returning
220  * the area for the rest of the earth.
221  * @param[out] perimeter the approximate perimeter of the polygon or length
222  * of the polyline (meters).
223  * @param[out] area the approximate area of the polygon
224  * (meters<sup>2</sup>); only set if polyline is false in the
225  * constructor.
226  * @return the number of points.
227  **********************************************************************/
228  unsigned TestEdge(real azi, real s, bool reverse, bool sign,
229  real& perimeter, real& area) const;
230 
231  /** \name Inspector functions
232  **********************************************************************/
233  ///@{
234  /**
235  * @return \e a the equatorial radius of the ellipsoid (meters). This is
236  * the value inherited from the Geodesic object used in the constructor.
237  **********************************************************************/
238 
239  Math::real EquatorialRadius() const { return _earth.EquatorialRadius(); }
240 
241  /**
242  * @return \e f the flattening of the ellipsoid. This is the value
243  * inherited from the Geodesic object used in the constructor.
244  **********************************************************************/
245  Math::real Flattening() const { return _earth.Flattening(); }
246 
247  /**
248  * Report the previous vertex added to the polygon or polyline.
249  *
250  * @param[out] lat the latitude of the point (degrees).
251  * @param[out] lon the longitude of the point (degrees).
252  *
253  * If no points have been added, then NaNs are returned. Otherwise, \e lon
254  * will be in the range [&minus;180&deg;, 180&deg;].
255  **********************************************************************/
256  void CurrentPoint(real& lat, real& lon) const
257  { lat = _lat1; lon = _lon1; }
258 
259  /**
260  * \deprecated An old name for EquatorialRadius().
261  **********************************************************************/
262  // GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
264  ///@}
265  };
266 
267  /**
268  * @relates PolygonAreaT
269  *
270  * Polygon areas using Geodesic. This should be used if the flattening is
271  * small.
272  **********************************************************************/
274 
275  /**
276  * @relates PolygonAreaT
277  *
278  * Polygon areas using GeodesicExact. (But note that the implementation of
279  * areas in GeodesicExact uses a high order series and this is only accurate
280  * for modest flattenings.)
281  **********************************************************************/
283 
284  /**
285  * @relates PolygonAreaT
286  *
287  * Polygon areas using Rhumb.
288  **********************************************************************/
290 
291 } // namespace GeographicLib
292 
293 #endif // GEOGRAPHICLIB_POLYGONAREA_HPP
static T AngNormalize(T x)
Definition: Math.hpp:383
Math::real Flattening() const
unsigned TestEdge(real azi, real s, bool reverse, bool sign, real &perimeter, real &area) const
Math::real EquatorialRadius() const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
PolygonAreaT< Rhumb > PolygonAreaRhumb
Accumulator & remainder(T y)
static T AngDiff(T x, T y, T &e)
Definition: Math.hpp:414
Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes.
void AddEdge(real azi, real s)
Definition: PolygonArea.cpp:38
PolygonAreaT(const GeodType &earth, bool polyline=false)
An accumulator for sums.
Definition: Accumulator.hpp:40
Header for GeographicLib::Geodesic class.
PolygonAreaT< GeodesicExact > PolygonAreaExact
Math::real MajorRadius() const
Header for GeographicLib::Accumulator class.
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
unsigned TestPoint(real lat, real lon, bool reverse, bool sign, real &perimeter, real &area) const
Definition: PolygonArea.cpp:81
static T NaN()
Definition: Math.cpp:389
unsigned Compute(bool reverse, bool sign, real &perimeter, real &area) const
Definition: PolygonArea.cpp:55
void AddPoint(real lat, real lon)
Definition: PolygonArea.cpp:17
void CurrentPoint(real &lat, real &lon) const
Header for GeographicLib::GeodesicExact class.
PolygonAreaT< Geodesic > PolygonArea
static T remainder(T x, T y)
Definition: Math.cpp:140