1 // Boost.Geometry (aka GGL, Generic Geometry Library) 2 3 // Copyright (c) 2016-2017 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 7 // Use, modification and distribution is subject to the Boost Software License, 8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 9 // http://www.boost.org/LICENSE_1_0.txt) 10 11 #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP 12 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP 13 14 15 #include <boost/geometry/core/srs.hpp> 16 17 #include <boost/geometry/formulas/area_formulas.hpp> 18 #include <boost/geometry/formulas/flattening.hpp> 19 20 #include <boost/geometry/strategies/geographic/parameters.hpp> 21 22 #include <boost/math/special_functions/atanh.hpp> 23 24 25 namespace boost { namespace geometry 26 { 27 28 namespace strategy { namespace area 29 { 30 31 /*! 32 \brief Geographic area calculation 33 \ingroup strategies 34 \details Geographic area calculation by trapezoidal rule plus integral 35 approximation that gives the ellipsoidal correction 36 \tparam PointOfSegment \tparam_segment_point 37 \tparam FormulaPolicy Formula used to calculate azimuths 38 \tparam SeriesOrder The order of approximation of the geodesic integral 39 \tparam Spheroid The spheroid model 40 \tparam CalculationType \tparam_calculation 41 \author See 42 - Danielsen JS, The area under the geodesic. Surv Rev 30(232): 61–66, 1989 43 - Charles F.F Karney, Algorithms for geodesics, 2011 https://arxiv.org/pdf/1109.4448.pdf 44 45 \qbk{ 46 [heading See also] 47 [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)] 48 } 49 */ 50 template 51 < 52 typename PointOfSegment, 53 typename FormulaPolicy = strategy::andoyer, 54 std::size_t SeriesOrder = strategy::default_order<FormulaPolicy>::value, 55 typename Spheroid = srs::spheroid<double>, 56 typename CalculationType = void 57 > 58 class geographic 59 { 60 // Switch between two kinds of approximation(series in eps and n v.s.series in k ^ 2 and e'^2) 61 static const bool ExpandEpsN = true; 62 // LongSegment Enables special handling of long segments 63 static const bool LongSegment = false; 64 65 //Select default types in case they are not set 66 67 typedef typename boost::mpl::if_c 68 < 69 boost::is_void<CalculationType>::type::value, 70 typename select_most_precise 71 < 72 typename coordinate_type<PointOfSegment>::type, 73 double 74 >::type, 75 CalculationType 76 >::type CT; 77 78 protected : 79 struct spheroid_constants 80 { 81 Spheroid m_spheroid; 82 CT const m_a2; // squared equatorial radius 83 CT const m_e2; // squared eccentricity 84 CT const m_ep2; // squared second eccentricity 85 CT const m_ep; // second eccentricity 86 CT const m_c2; // squared authalic radius 87 spheroid_constantsboost::geometry::strategy::area::geographic::spheroid_constants88 inline spheroid_constants(Spheroid const& spheroid) 89 : m_spheroid(spheroid) 90 , m_a2(math::sqr(get_radius<0>(spheroid))) 91 , m_e2(formula::flattening<CT>(spheroid) 92 * (CT(2.0) - CT(formula::flattening<CT>(spheroid)))) 93 , m_ep2(m_e2 / (CT(1.0) - m_e2)) 94 , m_ep(math::sqrt(m_ep2)) 95 , m_c2(authalic_radius(spheroid, m_a2, m_e2)) 96 {} 97 }; 98 authalic_radius(Spheroid const & sph,CT const & a2,CT const & e2)99 static inline CT authalic_radius(Spheroid const& sph, CT const& a2, CT const& e2) 100 { 101 CT const c0 = 0; 102 103 if (math::equals(e2, c0)) 104 { 105 return a2; 106 } 107 108 CT const sqrt_e2 = math::sqrt(e2); 109 CT const c2 = 2; 110 111 return (a2 / c2) + 112 ((math::sqr(get_radius<2>(sph)) * boost::math::atanh(sqrt_e2)) 113 / (c2 * sqrt_e2)); 114 } 115 116 struct area_sums 117 { 118 CT m_excess_sum; 119 CT m_correction_sum; 120 121 // Keep track if encircles some pole 122 std::size_t m_crosses_prime_meridian; 123 area_sumsboost::geometry::strategy::area::geographic::area_sums124 inline area_sums() 125 : m_excess_sum(0) 126 , m_correction_sum(0) 127 , m_crosses_prime_meridian(0) 128 {} areaboost::geometry::strategy::area::geographic::area_sums129 inline CT area(spheroid_constants spheroid_const) const 130 { 131 CT result; 132 133 CT sum = spheroid_const.m_c2 * m_excess_sum 134 + spheroid_const.m_e2 * spheroid_const.m_a2 * m_correction_sum; 135 136 // If encircles some pole 137 if (m_crosses_prime_meridian % 2 == 1) 138 { 139 std::size_t times_crosses_prime_meridian 140 = 1 + (m_crosses_prime_meridian / 2); 141 142 result = CT(2.0) 143 * geometry::math::pi<CT>() 144 * spheroid_const.m_c2 145 * CT(times_crosses_prime_meridian) 146 - geometry::math::abs(sum); 147 148 if (geometry::math::sign<CT>(sum) == 1) 149 { 150 result = - result; 151 } 152 153 } 154 else 155 { 156 result = sum; 157 } 158 159 return result; 160 } 161 }; 162 163 public : 164 typedef CT return_type; 165 typedef PointOfSegment segment_point_type; 166 typedef area_sums state_type; 167 geographic(Spheroid const & spheroid=Spheroid ())168 explicit inline geographic(Spheroid const& spheroid = Spheroid()) 169 : m_spheroid_constants(spheroid) 170 {} 171 apply(PointOfSegment const & p1,PointOfSegment const & p2,area_sums & state) const172 inline void apply(PointOfSegment const& p1, 173 PointOfSegment const& p2, 174 area_sums& state) const 175 { 176 177 if (! geometry::math::equals(get<0>(p1), get<0>(p2))) 178 { 179 180 typedef geometry::formula::area_formulas 181 < 182 CT, SeriesOrder, ExpandEpsN 183 > area_formulas; 184 185 typename area_formulas::return_type_ellipsoidal result = 186 area_formulas::template ellipsoidal<FormulaPolicy::template inverse> 187 (p1, p2, m_spheroid_constants); 188 189 state.m_excess_sum += result.spherical_term; 190 state.m_correction_sum += result.ellipsoidal_term; 191 192 // Keep track whenever a segment crosses the prime meridian 193 geometry::formula::area_formulas<CT> 194 ::crosses_prime_meridian(p1, p2, state); 195 } 196 } 197 result(area_sums const & state) const198 inline return_type result(area_sums const& state) const 199 { 200 return state.area(m_spheroid_constants); 201 } 202 203 private: 204 spheroid_constants m_spheroid_constants; 205 206 }; 207 208 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 209 210 namespace services 211 { 212 213 214 template <typename Point> 215 struct default_strategy<geographic_tag, Point> 216 { 217 typedef strategy::area::geographic<Point> type; 218 }; 219 220 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 221 222 } 223 224 }} // namespace strategy::area 225 226 227 228 229 }} // namespace boost::geometry 230 231 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP 232