strategy/geographic/area.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2016-2020 Oracle and/or its affiliates.
0006 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 #ifndef BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_HPP
0014 #define BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_HPP
0015
0016
0017 #include <type_traits>
0018
0019 #include <boost/geometry/srs/spheroid.hpp>
0020
0021 #include <boost/geometry/formulas/area_formulas.hpp>
0022 #include <boost/geometry/formulas/authalic_radius_sqr.hpp>
0023 #include <boost/geometry/formulas/eccentricity_sqr.hpp>
0024
0025 #include <boost/geometry/strategy/area.hpp>
0026 #include <boost/geometry/strategies/geographic/parameters.hpp>
0027
0028
0029 namespace boost { namespace geometry
0030 {
0031
0032 namespace strategy { namespace area
0033 {
0034
0035 /*!
0036 \brief Geographic area calculation
0037 \ingroup strategies
0038 \details Geographic area calculation by trapezoidal rule plus integral
0039          approximation that gives the ellipsoidal correction
0040 \tparam FormulaPolicy Formula used to calculate azimuths
0041 \tparam SeriesOrder The order of approximation of the geodesic integral
0042 \tparam Spheroid The spheroid model
0043 \tparam CalculationType \tparam_calculation
0044 \author See
0045 - Danielsen JS, The area under the geodesic. Surv Rev 30(232): 61–66, 1989
0046 - Charles F.F Karney, Algorithms for geodesics, 2011 https://arxiv.org/pdf/1109.4448.pdf
0047
0048 \qbk{
0049 [heading See also]
0050 \* [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)]
0051 \* [link geometry.reference.srs.srs_spheroid srs::spheroid]
0052 }
0053 */
0054 template
0055 <
0056     typename FormulaPolicy = strategy::andoyer,
0057     std::size_t SeriesOrder = strategy::default_order<FormulaPolicy>::value,
0058     typename Spheroid = srs::spheroid<double>,
0059     typename CalculationType = void
0060 >
0061 class geographic
0062 {
0063     // Switch between two kinds of approximation(series in eps and n v.s.series in k ^ 2 and e'^2)
0064     static const bool ExpandEpsN = true;
0065     // LongSegment Enables special handling of long segments
0066     static const bool LongSegment = false;
0067     // Area formula is implemented for a maximum series order 5
0068     static constexpr auto SeriesOrderNorm = SeriesOrder > 5 ? 5 : SeriesOrder;
0069
0070     //Select default types in case they are not set
0071
0072 public:
0073     template <typename Geometry>
0074     struct result_type
0075         : strategy::area::detail::result_type
0076             <
0077                 Geometry,
0078                 CalculationType
0079             >
0080     {};
0081
0082 protected :
0083     struct spheroid_constants
0084     {
0085         typedef std::conditional_t
0086             <
0087                 std::is_void<CalculationType>::value,
0088                 typename geometry::radius_type<Spheroid>::type,
0089                 CalculationType
0090             > calc_t;
0091
0092         Spheroid m_spheroid;
0093         calc_t const m_a2;  // squared equatorial radius
0094         calc_t const m_e2;  // squared eccentricity
0095         calc_t const m_ep2; // squared second eccentricity
0096         calc_t const m_ep;  // second eccentricity
0097         calc_t const m_c2;  // squared authalic radius
0098         calc_t const m_f;   // the flattening
0099         calc_t m_coeffs_var[((SeriesOrderNorm+2)*(SeriesOrderNorm+1))/2];
0100
0101         inline spheroid_constants(Spheroid const& spheroid)
0102             : m_spheroid(spheroid)
0103             , m_a2(math::sqr(get_radius<0>(spheroid)))
0104             , m_e2(formula::eccentricity_sqr<calc_t>(spheroid))
0105             , m_ep2(m_e2 / (calc_t(1.0) - m_e2))
0106             , m_ep(math::sqrt(m_ep2))
0107             , m_c2(formula_dispatch::authalic_radius_sqr
0108                     <
0109                         calc_t, Spheroid, srs_spheroid_tag
0110                     >::apply(m_a2, m_e2))
0111             , m_f(formula::flattening<calc_t>(spheroid))
0112         {
0113             typedef geometry::formula::area_formulas
0114                 <
0115                     calc_t, SeriesOrderNorm, ExpandEpsN
0116                 > area_formulas;
0117
0118             calc_t const n = m_f / (calc_t(2) - m_f);
0119
0120             // Generate and evaluate the polynomials on n
0121             // to get the series coefficients (that depend on eps)
0122             area_formulas::evaluate_coeffs_n(n, m_coeffs_var);
0123         }
0124     };
0125
0126 public:
0127     template <typename Geometry>
0128     class state
0129     {
0130         friend class geographic;
0131
0132         typedef typename result_type<Geometry>::type return_type;
0133
0134     public:
0135         inline state()
0136             : m_excess_sum(0)
0137             , m_correction_sum(0)
0138             , m_crosses_prime_meridian(0)
0139         {}
0140
0141     private:
0142         inline return_type area(spheroid_constants const& spheroid_const) const
0143         {
0144             return_type result;
0145
0146             return_type const spherical_term = spheroid_const.m_c2 * m_excess_sum;
0147             return_type const ellipsoidal_term = spheroid_const.m_e2
0148                 * spheroid_const.m_a2 * m_correction_sum;
0149
0150             // ignore ellipsoidal term if is large (probably from an azimuth
0151             // inaccuracy)
0152             return_type sum = math::abs(ellipsoidal_term/spherical_term) > 0.01
0153                 ? spherical_term : spherical_term + ellipsoidal_term;
0154
0155             // If encircles some pole
0156             if (m_crosses_prime_meridian % 2 == 1)
0157             {
0158                 std::size_t times_crosses_prime_meridian
0159                         = 1 + (m_crosses_prime_meridian / 2);
0160
0161                 result = return_type(2.0)
0162                          * geometry::math::pi<return_type>()
0163                          * spheroid_const.m_c2
0164                          * return_type(times_crosses_prime_meridian)
0165                          - geometry::math::abs(sum);
0166
0167                 if (geometry::math::sign<return_type>(sum) == 1)
0168                 {
0169                     result = - result;
0170                 }
0171
0172             }
0173             else
0174             {
0175                 result = sum;
0176             }
0177
0178             return result;
0179         }
0180
0181         return_type m_excess_sum;
0182         return_type m_correction_sum;
0183
0184         // Keep track if encircles some pole
0185         std::size_t m_crosses_prime_meridian;
0186     };
0187
0188 public :
0189     explicit inline geographic(Spheroid const& spheroid = Spheroid())
0190         : m_spheroid_constants(spheroid)
0191     {}
0192
0193     template <typename PointOfSegment, typename Geometry>
0194     inline void apply(PointOfSegment const& p1,
0195                       PointOfSegment const& p2,
0196                       state<Geometry>& st) const
0197     {
0198         using CT = typename result_type<Geometry>::type;
0199
0200         // if the segment in not on a meridian
0201         if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
0202         {
0203             typedef geometry::formula::area_formulas
0204                 <
0205                     CT, SeriesOrderNorm, ExpandEpsN
0206                 > area_formulas;
0207
0208             // Keep track whenever a segment crosses the prime meridian
0209             if (area_formulas::crosses_prime_meridian(p1, p2))
0210             {
0211                 st.m_crosses_prime_meridian++;
0212             }
0213
0214             // if the segment in not on equator
0215             if (! (geometry::math::equals(get<1>(p1), 0)
0216                 && geometry::math::equals(get<1>(p2), 0)))
0217             {
0218                 auto result = area_formulas::template ellipsoidal
0219                     <
0220                         FormulaPolicy::template inverse
0221                     >(p1, p2, m_spheroid_constants);
0222
0223                 st.m_excess_sum += result.spherical_term;
0224                 st.m_correction_sum += result.ellipsoidal_term;
0225             }
0226         }
0227     }
0228
0229     template <typename Geometry>
0230     inline typename result_type<Geometry>::type
0231         result(state<Geometry> const& st) const
0232     {
0233         return st.area(m_spheroid_constants);
0234     }
0235
0236     Spheroid model() const
0237     {
0238         return m_spheroid_constants.m_spheroid;
0239     }
0240
0241 private:
0242     spheroid_constants m_spheroid_constants;
0243
0244 };
0245
0246 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0247
0248 namespace services
0249 {
0250
0251
0252 template <>
0253 struct default_strategy<geographic_tag>
0254 {
0255     typedef strategy::area::geographic<> type;
0256 };
0257
0258 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0259
0260 }
0261
0262 }} // namespace strategy::area
0263
0264
0265
0266
0267 }} // namespace boost::geometry
0268
0269 #endif // BOOST_GEOMETRY_STRATEGY_GEOGRAPHIC_AREA_HPP