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