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