GeographicLib  1.49
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-2016) <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  * This class lets you add vertices and edges one at a time to the polygon.
36  * The sequence must start with a vertex and thereafter vertices and edges
37  * can be added in any order. Any vertex after the first creates a new edge
38  * which is the \e shortest geodesic from the previous vertex. In some
39  * cases there may be two or many such shortest geodesics and the area is
40  * then not uniquely defined. In this case, either add an intermediate
41  * vertex or add the edge \e as an edge (by defining its direction and
42  * length).
43  *
44  * The area and perimeter are accumulated at two times the standard floating
45  * point precision to guard against the loss of accuracy with many-sided
46  * polygons. At any point you can ask for the perimeter and area so far.
47  * There's an option to treat the points as defining a polyline instead of a
48  * polygon; in that case, only the perimeter is computed.
49  *
50  * This is a templated class to allow it to be used with Geodesic,
51  * GeodesicExact, and Rhumb. GeographicLib::PolygonArea,
52  * GeographicLib::PolygonAreaExact, and GeographicLib::PolygonAreaRhumb are
53  * typedefs for these cases.
54  *
55  * @tparam GeodType the geodesic class to use.
56  *
57  * Example of use:
58  * \include example-PolygonArea.cpp
59  *
60  * <a href="Planimeter.1.html">Planimeter</a> is a command-line utility
61  * providing access to the functionality of PolygonAreaT.
62  **********************************************************************/
63 
64  template <class GeodType = Geodesic>
65  class PolygonAreaT {
66  private:
67  typedef Math::real real;
68  GeodType _earth;
69  real _area0; // Full ellipsoid area
70  bool _polyline; // Assume polyline (don't close and skip area)
71  unsigned _mask;
72  unsigned _num;
73  int _crossings;
74  Accumulator<> _areasum, _perimetersum;
75  real _lat0, _lon0, _lat1, _lon1;
76  static int transit(real lon1, real lon2) {
77  // Return 1 or -1 if crossing prime meridian in east or west direction.
78  // Otherwise return zero.
79  // Compute lon12 the same way as Geodesic::Inverse.
80  lon1 = Math::AngNormalize(lon1);
81  lon2 = Math::AngNormalize(lon2);
82  real lon12 = Math::AngDiff(lon1, lon2);
83  // Treat 0 as negative in these tests. This balances +/- 180 being
84  // treated as positive, i.e., +180.
85  int cross =
86  lon1 <= 0 && lon2 > 0 && lon12 > 0 ? 1 :
87  (lon2 <= 0 && lon1 > 0 && lon12 < 0 ? -1 : 0);
88  return cross;
89  }
90  // an alternate version of transit to deal with longitudes in the direct
91  // problem.
92  static int transitdirect(real lon1, real lon2) {
93  // We want to compute exactly
94  // int(floor(lon2 / 360)) - int(floor(lon1 / 360))
95  // Since we only need the parity of the result we can use std::remquo;
96  // but this is buggy with g++ 4.8.3 (glibc version < 2.22), see
97  // https://sourceware.org/bugzilla/show_bug.cgi?id=17569
98  // and requires C++11. So instead we do
99 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 4
100  using std::remainder;
101  lon1 = remainder(lon1, real(720)); lon2 = remainder(lon2, real(720));
102  return ( (lon2 >= 0 && lon2 < 360 ? 0 : 1) -
103  (lon1 >= 0 && lon1 < 360 ? 0 : 1) );
104 #else
105  using std::fmod;
106  lon1 = fmod(lon1, real(720)); lon2 = fmod(lon2, real(720));
107  return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
108  ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
109 #endif
110  }
111  public:
112 
113  /**
114  * Constructor for PolygonAreaT.
115  *
116  * @param[in] earth the Geodesic object to use for geodesic calculations.
117  * @param[in] polyline if true that treat the points as defining a polyline
118  * instead of a polygon (default = false).
119  **********************************************************************/
120  PolygonAreaT(const GeodType& earth, bool polyline = false)
121  : _earth(earth)
122  , _area0(_earth.EllipsoidArea())
123  , _polyline(polyline)
124  , _mask(GeodType::LATITUDE | GeodType::LONGITUDE | GeodType::DISTANCE |
125  (_polyline ? GeodType::NONE :
126  GeodType::AREA | GeodType::LONG_UNROLL))
127  { Clear(); }
128 
129  /**
130  * Clear PolygonAreaT, allowing a new polygon to be started.
131  **********************************************************************/
132  void Clear() {
133  _num = 0;
134  _crossings = 0;
135  _areasum = 0;
136  _perimetersum = 0;
137  _lat0 = _lon0 = _lat1 = _lon1 = Math::NaN();
138  }
139 
140  /**
141  * Add a point to the polygon or polyline.
142  *
143  * @param[in] lat the latitude of the point (degrees).
144  * @param[in] lon the longitude of the point (degrees).
145  *
146  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
147  **********************************************************************/
148  void AddPoint(real lat, real lon);
149 
150  /**
151  * Add an edge to the polygon or polyline.
152  *
153  * @param[in] azi azimuth at current point (degrees).
154  * @param[in] s distance from current point to next point (meters).
155  *
156  * This does nothing if no points have been added yet. Use
157  * PolygonAreaT::CurrentPoint to determine the position of the new vertex.
158  **********************************************************************/
159  void AddEdge(real azi, real s);
160 
161  /**
162  * Return the results so far.
163  *
164  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
165  * traversal counts as a positive area.
166  * @param[in] sign if true then return a signed result for the area if
167  * the polygon is traversed in the "wrong" direction instead of returning
168  * the area for the rest of the earth.
169  * @param[out] perimeter the perimeter of the polygon or length of the
170  * polyline (meters).
171  * @param[out] area the area of the polygon (meters<sup>2</sup>); only set
172  * if \e polyline is false in the constructor.
173  * @return the number of points.
174  *
175  * More points can be added to the polygon after this call.
176  **********************************************************************/
177  unsigned Compute(bool reverse, bool sign,
178  real& perimeter, real& area) const;
179 
180  /**
181  * Return the results assuming a tentative final test point is added;
182  * however, the data for the test point is not saved. This lets you report
183  * a running result for the perimeter and area as the user moves the mouse
184  * cursor. Ordinary floating point arithmetic is used to accumulate the
185  * data for the test point; thus the area and perimeter returned are less
186  * accurate than if PolygonAreaT::AddPoint and PolygonAreaT::Compute are
187  * used.
188  *
189  * @param[in] lat the latitude of the test point (degrees).
190  * @param[in] lon the longitude of the test point (degrees).
191  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
192  * traversal counts as a positive area.
193  * @param[in] sign if true then return a signed result for the area if
194  * the polygon is traversed in the "wrong" direction instead of returning
195  * the area for the rest of the earth.
196  * @param[out] perimeter the approximate perimeter of the polygon or length
197  * of the polyline (meters).
198  * @param[out] area the approximate area of the polygon
199  * (meters<sup>2</sup>); only set if polyline is false in the
200  * constructor.
201  * @return the number of points.
202  *
203  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
204  **********************************************************************/
205  unsigned TestPoint(real lat, real lon, bool reverse, bool sign,
206  real& perimeter, real& area) const;
207 
208  /**
209  * Return the results assuming a tentative final test point is added via an
210  * azimuth and distance; however, the data for the test point is not saved.
211  * This lets you report a running result for the perimeter and area as the
212  * user moves the mouse cursor. Ordinary floating point arithmetic is used
213  * to accumulate the data for the test point; thus the area and perimeter
214  * returned are less accurate than if PolygonAreaT::AddEdge and
215  * PolygonAreaT::Compute are used.
216  *
217  * @param[in] azi azimuth at current point (degrees).
218  * @param[in] s distance from current point to final test point (meters).
219  * @param[in] reverse if true then clockwise (instead of counter-clockwise)
220  * traversal counts as a positive area.
221  * @param[in] sign if true then return a signed result for the area if
222  * the polygon is traversed in the "wrong" direction instead of returning
223  * the area for the rest of the earth.
224  * @param[out] perimeter the approximate perimeter of the polygon or length
225  * of the polyline (meters).
226  * @param[out] area the approximate area of the polygon
227  * (meters<sup>2</sup>); only set if polyline is false in the
228  * constructor.
229  * @return the number of points.
230  **********************************************************************/
231  unsigned TestEdge(real azi, real s, bool reverse, bool sign,
232  real& perimeter, real& area) const;
233 
234  /** \name Inspector functions
235  **********************************************************************/
236  ///@{
237  /**
238  * @return \e a the equatorial radius of the ellipsoid (meters). This is
239  * the value inherited from the Geodesic object used in the constructor.
240  **********************************************************************/
241 
242  Math::real MajorRadius() const { return _earth.MajorRadius(); }
243 
244  /**
245  * @return \e f the flattening of the ellipsoid. This is the value
246  * inherited from the Geodesic object used in the constructor.
247  **********************************************************************/
248  Math::real Flattening() const { return _earth.Flattening(); }
249 
250  /**
251  * Report the previous vertex added to the polygon or polyline.
252  *
253  * @param[out] lat the latitude of the point (degrees).
254  * @param[out] lon the longitude of the point (degrees).
255  *
256  * If no points have been added, then NaNs are returned. Otherwise, \e lon
257  * will be in the range [&minus;180&deg;, 180&deg;].
258  **********************************************************************/
259  void CurrentPoint(real& lat, real& lon) const
260  { lat = _lat1; lon = _lon1; }
261  ///@}
262  };
263 
264  /**
265  * @relates PolygonAreaT
266  *
267  * Polygon areas using Geodesic. This should be used if the flattening is
268  * small.
269  **********************************************************************/
271 
272  /**
273  * @relates PolygonAreaT
274  *
275  * Polygon areas using GeodesicExact. (But note that the implementation of
276  * areas in GeodesicExact uses a high order series and this is only accurate
277  * for modest flattenings.)
278  **********************************************************************/
280 
281  /**
282  * @relates PolygonAreaT
283  *
284  * Polygon areas using Rhumb.
285  **********************************************************************/
287 
288 } // namespace GeographicLib
289 
290 #endif // GEOGRAPHICLIB_POLYGONAREA_HPP
static T AngNormalize(T x)
Definition: Math.hpp:440
static T NaN()
Definition: Math.hpp:830
Math::real Flattening() const
unsigned TestEdge(real azi, real s, bool reverse, bool sign, real &perimeter, real &area) const
PolygonAreaT< Rhumb > PolygonAreaRhumb
static T AngDiff(T x, T y, T &e)
Definition: Math.hpp:486
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:98
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