area.hpp 6.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
  3. // Copyright (c) 2016-2018 Oracle and/or its affiliates.
  4. // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
  5. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  6. // Use, modification and distribution is subject to the Boost Software License,
  7. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8. // http://www.boost.org/LICENSE_1_0.txt)
  9. #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP
  10. #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP
  11. #include <boost/geometry/srs/spheroid.hpp>
  12. #include <boost/geometry/formulas/area_formulas.hpp>
  13. #include <boost/geometry/formulas/authalic_radius_sqr.hpp>
  14. #include <boost/geometry/formulas/eccentricity_sqr.hpp>
  15. #include <boost/geometry/strategies/geographic/parameters.hpp>
  16. namespace boost { namespace geometry
  17. {
  18. namespace strategy { namespace area
  19. {
  20. /*!
  21. \brief Geographic area calculation
  22. \ingroup strategies
  23. \details Geographic area calculation by trapezoidal rule plus integral
  24. approximation that gives the ellipsoidal correction
  25. \tparam FormulaPolicy Formula used to calculate azimuths
  26. \tparam SeriesOrder The order of approximation of the geodesic integral
  27. \tparam Spheroid The spheroid model
  28. \tparam CalculationType \tparam_calculation
  29. \author See
  30. - Danielsen JS, The area under the geodesic. Surv Rev 30(232): 61–66, 1989
  31. - Charles F.F Karney, Algorithms for geodesics, 2011 https://arxiv.org/pdf/1109.4448.pdf
  32. \qbk{
  33. [heading See also]
  34. \* [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)]
  35. \* [link geometry.reference.srs.srs_spheroid srs::spheroid]
  36. }
  37. */
  38. template
  39. <
  40. typename FormulaPolicy = strategy::andoyer,
  41. std::size_t SeriesOrder = strategy::default_order<FormulaPolicy>::value,
  42. typename Spheroid = srs::spheroid<double>,
  43. typename CalculationType = void
  44. >
  45. class geographic
  46. {
  47. // Switch between two kinds of approximation(series in eps and n v.s.series in k ^ 2 and e'^2)
  48. static const bool ExpandEpsN = true;
  49. // LongSegment Enables special handling of long segments
  50. static const bool LongSegment = false;
  51. //Select default types in case they are not set
  52. public:
  53. template <typename Geometry>
  54. struct result_type
  55. : strategy::area::detail::result_type
  56. <
  57. Geometry,
  58. CalculationType
  59. >
  60. {};
  61. protected :
  62. struct spheroid_constants
  63. {
  64. typedef typename boost::mpl::if_c
  65. <
  66. boost::is_void<CalculationType>::value,
  67. typename geometry::radius_type<Spheroid>::type,
  68. CalculationType
  69. >::type calc_t;
  70. Spheroid m_spheroid;
  71. calc_t const m_a2; // squared equatorial radius
  72. calc_t const m_e2; // squared eccentricity
  73. calc_t const m_ep2; // squared second eccentricity
  74. calc_t const m_ep; // second eccentricity
  75. calc_t const m_c2; // squared authalic radius
  76. inline spheroid_constants(Spheroid const& spheroid)
  77. : m_spheroid(spheroid)
  78. , m_a2(math::sqr(get_radius<0>(spheroid)))
  79. , m_e2(formula::eccentricity_sqr<calc_t>(spheroid))
  80. , m_ep2(m_e2 / (calc_t(1.0) - m_e2))
  81. , m_ep(math::sqrt(m_ep2))
  82. , m_c2(formula_dispatch::authalic_radius_sqr
  83. <
  84. calc_t, Spheroid, srs_spheroid_tag
  85. >::apply(m_a2, m_e2))
  86. {}
  87. };
  88. public:
  89. template <typename Geometry>
  90. class state
  91. {
  92. friend class geographic;
  93. typedef typename result_type<Geometry>::type return_type;
  94. public:
  95. inline state()
  96. : m_excess_sum(0)
  97. , m_correction_sum(0)
  98. , m_crosses_prime_meridian(0)
  99. {}
  100. private:
  101. inline return_type area(spheroid_constants const& spheroid_const) const
  102. {
  103. return_type result;
  104. return_type sum = spheroid_const.m_c2 * m_excess_sum
  105. + spheroid_const.m_e2 * spheroid_const.m_a2 * m_correction_sum;
  106. // If encircles some pole
  107. if (m_crosses_prime_meridian % 2 == 1)
  108. {
  109. std::size_t times_crosses_prime_meridian
  110. = 1 + (m_crosses_prime_meridian / 2);
  111. result = return_type(2.0)
  112. * geometry::math::pi<return_type>()
  113. * spheroid_const.m_c2
  114. * return_type(times_crosses_prime_meridian)
  115. - geometry::math::abs(sum);
  116. if (geometry::math::sign<return_type>(sum) == 1)
  117. {
  118. result = - result;
  119. }
  120. }
  121. else
  122. {
  123. result = sum;
  124. }
  125. return result;
  126. }
  127. return_type m_excess_sum;
  128. return_type m_correction_sum;
  129. // Keep track if encircles some pole
  130. std::size_t m_crosses_prime_meridian;
  131. };
  132. public :
  133. explicit inline geographic(Spheroid const& spheroid = Spheroid())
  134. : m_spheroid_constants(spheroid)
  135. {}
  136. template <typename PointOfSegment, typename Geometry>
  137. inline void apply(PointOfSegment const& p1,
  138. PointOfSegment const& p2,
  139. state<Geometry>& st) const
  140. {
  141. if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
  142. {
  143. typedef geometry::formula::area_formulas
  144. <
  145. typename result_type<Geometry>::type,
  146. SeriesOrder, ExpandEpsN
  147. > area_formulas;
  148. typename area_formulas::return_type_ellipsoidal result =
  149. area_formulas::template ellipsoidal<FormulaPolicy::template inverse>
  150. (p1, p2, m_spheroid_constants);
  151. st.m_excess_sum += result.spherical_term;
  152. st.m_correction_sum += result.ellipsoidal_term;
  153. // Keep track whenever a segment crosses the prime meridian
  154. if (area_formulas::crosses_prime_meridian(p1, p2))
  155. {
  156. st.m_crosses_prime_meridian++;
  157. }
  158. }
  159. }
  160. template <typename Geometry>
  161. inline typename result_type<Geometry>::type
  162. result(state<Geometry> const& st) const
  163. {
  164. return st.area(m_spheroid_constants);
  165. }
  166. private:
  167. spheroid_constants m_spheroid_constants;
  168. };
  169. #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  170. namespace services
  171. {
  172. template <>
  173. struct default_strategy<geographic_tag>
  174. {
  175. typedef strategy::area::geographic<> type;
  176. };
  177. #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  178. }
  179. }} // namespace strategy::area
  180. }} // namespace boost::geometry
  181. #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP