]>
Commit | Line | Data |
---|---|---|
92f5a8d4 | 1 | // Boost.Geometry - gis-projections (based on PROJ4) |
11fdf7f2 TL |
2 | |
3 | // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. | |
4 | ||
92f5a8d4 TL |
5 | // This file was modified by Oracle on 2017, 2018, 2019. |
6 | // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates. | |
11fdf7f2 TL |
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 | // This file is converted from PROJ4, http://trac.osgeo.org/proj | |
14 | // PROJ4 is originally written by Gerald Evenden (then of the USGS) | |
15 | // PROJ4 is maintained by Frank Warmerdam | |
16 | // PROJ4 is converted to Boost.Geometry by Barend Gehrels | |
17 | ||
92f5a8d4 | 18 | // Last updated version of proj: 5.0.0 |
11fdf7f2 TL |
19 | |
20 | // Original copyright notice: | |
21 | ||
22 | // Purpose: Implementation of the aitoff (Aitoff) and wintri (Winkel Tripel) | |
92f5a8d4 TL |
23 | // projections. |
24 | // Author: Gerald Evenden (1995) | |
25 | // Drazen Tutic, Lovro Gradiser (2015) - add inverse | |
26 | // Thomas Knudsen (2016) - revise/add regression tests | |
11fdf7f2 TL |
27 | // Copyright (c) 1995, Gerald Evenden |
28 | ||
29 | // Permission is hereby granted, free of charge, to any person obtaining a | |
30 | // copy of this software and associated documentation files (the "Software"), | |
31 | // to deal in the Software without restriction, including without limitation | |
32 | // the rights to use, copy, modify, merge, publish, distribute, sublicense, | |
33 | // and/or sell copies of the Software, and to permit persons to whom the | |
34 | // Software is furnished to do so, subject to the following conditions: | |
35 | ||
36 | // The above copyright notice and this permission notice shall be included | |
37 | // in all copies or substantial portions of the Software. | |
38 | ||
39 | // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS | |
40 | // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | |
41 | // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL | |
42 | // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | |
43 | // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING | |
44 | // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER | |
45 | // DEALINGS IN THE SOFTWARE. | |
46 | ||
92f5a8d4 TL |
47 | #ifndef BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP |
48 | #define BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP | |
49 | ||
11fdf7f2 | 50 | #include <boost/core/ignore_unused.hpp> |
11fdf7f2 TL |
51 | |
52 | #include <boost/geometry/srs/projections/impl/base_static.hpp> | |
53 | #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> | |
11fdf7f2 | 54 | #include <boost/geometry/srs/projections/impl/factory_entry.hpp> |
92f5a8d4 TL |
55 | #include <boost/geometry/srs/projections/impl/pj_param.hpp> |
56 | #include <boost/geometry/srs/projections/impl/projects.hpp> | |
11fdf7f2 | 57 | |
92f5a8d4 | 58 | #include <boost/geometry/util/math.hpp> |
11fdf7f2 | 59 | |
92f5a8d4 | 60 | namespace boost { namespace geometry |
11fdf7f2 | 61 | { |
11fdf7f2 TL |
62 | |
63 | namespace projections | |
64 | { | |
65 | #ifndef DOXYGEN_NO_DETAIL | |
66 | namespace detail { namespace aitoff | |
67 | { | |
92f5a8d4 TL |
68 | enum mode_type { |
69 | mode_aitoff = 0, | |
70 | mode_winkel_tripel = 1 | |
71 | }; | |
72 | ||
11fdf7f2 TL |
73 | template <typename T> |
74 | struct par_aitoff | |
75 | { | |
76 | T cosphi1; | |
92f5a8d4 | 77 | mode_type mode; |
11fdf7f2 TL |
78 | }; |
79 | ||
92f5a8d4 TL |
80 | template <typename T, typename Parameters> |
81 | struct base_aitoff_spheroid | |
11fdf7f2 | 82 | { |
92f5a8d4 | 83 | par_aitoff<T> m_proj_parm; |
11fdf7f2 TL |
84 | |
85 | // FORWARD(s_forward) spheroid | |
86 | // Project coordinates from geographic (lon, lat) to cartesian (x, y) | |
92f5a8d4 | 87 | inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const |
11fdf7f2 | 88 | { |
92f5a8d4 | 89 | T c, d; |
11fdf7f2 TL |
90 | |
91 | if((d = acos(cos(lp_lat) * cos(c = 0.5 * lp_lon)))) {/* basic Aitoff */ | |
92 | xy_x = 2. * d * cos(lp_lat) * sin(c) * (xy_y = 1. / sin(d)); | |
93 | xy_y *= d * sin(lp_lat); | |
94 | } else | |
95 | xy_x = xy_y = 0.; | |
92f5a8d4 | 96 | if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */ |
11fdf7f2 TL |
97 | xy_x = (xy_x + lp_lon * this->m_proj_parm.cosphi1) * 0.5; |
98 | xy_y = (xy_y + lp_lat) * 0.5; | |
99 | } | |
100 | } | |
101 | /*********************************************************************************** | |
102 | * | |
103 | * Inverse functions added by Drazen Tutic and Lovro Gradiser based on paper: | |
104 | * | |
105 | * I.Özbug Biklirici and Cengizhan Ipbüker. A General Algorithm for the Inverse | |
106 | * Transformation of Map Projections Using Jacobian Matrices. In Proceedings of the | |
107 | * Third International Symposium Mathematical & Computational Applications, | |
108 | * pages 175{182, Turkey, September 2002. | |
109 | * | |
92f5a8d4 | 110 | * Expected accuracy is defined by epsilon = 1e-12. Should be appropriate for |
11fdf7f2 TL |
111 | * most applications of Aitoff and Winkel Tripel projections. |
112 | * | |
113 | * Longitudes of 180W and 180E can be mixed in solution obtained. | |
114 | * | |
115 | * Inverse for Aitoff projection in poles is undefined, longitude value of 0 is assumed. | |
116 | * | |
117 | * Contact : dtutic@geof.hr | |
118 | * Date: 2015-02-16 | |
119 | * | |
120 | ************************************************************************************/ | |
121 | ||
122 | // INVERSE(s_inverse) sphere | |
123 | // Project coordinates from cartesian (x, y) to geographic (lon, lat) | |
92f5a8d4 | 124 | inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const |
11fdf7f2 | 125 | { |
92f5a8d4 TL |
126 | static const T pi = detail::pi<T>(); |
127 | static const T two_pi = detail::two_pi<T>(); | |
128 | static const T epsilon = 1e-12; | |
11fdf7f2 | 129 | |
92f5a8d4 TL |
130 | int iter, max_iter = 10, round = 0, max_round = 20; |
131 | T D, C, f1, f2, f1p, f1l, f2p, f2l, dp, dl, sl, sp, cp, cl, x, y; | |
11fdf7f2 | 132 | |
92f5a8d4 TL |
133 | if ((fabs(xy_x) < epsilon) && (fabs(xy_y) < epsilon )) { |
134 | lp_lat = 0.; lp_lon = 0.; | |
135 | return; | |
136 | } | |
11fdf7f2 TL |
137 | |
138 | /* intial values for Newton-Raphson method */ | |
139 | lp_lat = xy_y; lp_lon = xy_x; | |
140 | do { | |
141 | iter = 0; | |
142 | do { | |
143 | sl = sin(lp_lon * 0.5); cl = cos(lp_lon * 0.5); | |
144 | sp = sin(lp_lat); cp = cos(lp_lat); | |
145 | D = cp * cl; | |
92f5a8d4 TL |
146 | C = 1. - D * D; |
147 | D = acos(D) / math::pow(C, T(1.5)); | |
148 | f1 = 2. * D * C * cp * sl; | |
149 | f2 = D * C * sp; | |
150 | f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl); | |
151 | f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp; | |
152 | f2p = sp * sp * cl / C + D * sl * sl * cp; | |
153 | f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl); | |
154 | if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */ | |
11fdf7f2 TL |
155 | f1 = 0.5 * (f1 + lp_lon * this->m_proj_parm.cosphi1); |
156 | f2 = 0.5 * (f2 + lp_lat); | |
157 | f1p *= 0.5; | |
158 | f1l = 0.5 * (f1l + this->m_proj_parm.cosphi1); | |
159 | f2p = 0.5 * (f2p + 1.); | |
160 | f2l *= 0.5; | |
161 | } | |
162 | f1 -= xy_x; f2 -= xy_y; | |
163 | dl = (f2 * f1p - f1 * f2p) / (dp = f1p * f2l - f2p * f1l); | |
164 | dp = (f1 * f2l - f2 * f1l) / dp; | |
92f5a8d4 | 165 | dl = fmod(dl, pi); /* set to interval [-M_PI, M_PI] */ |
11fdf7f2 | 166 | lp_lat -= dp; lp_lon -= dl; |
92f5a8d4 TL |
167 | } while ((fabs(dp) > epsilon || fabs(dl) > epsilon) && (iter++ < max_iter)); |
168 | if (lp_lat > two_pi) lp_lat -= 2.*(lp_lat-two_pi); /* correct if symmetrical solution for Aitoff */ | |
169 | if (lp_lat < -two_pi) lp_lat -= 2.*(lp_lat+two_pi); /* correct if symmetrical solution for Aitoff */ | |
170 | if ((fabs(fabs(lp_lat) - two_pi) < epsilon) && (!this->m_proj_parm.mode)) lp_lon = 0.; /* if pole in Aitoff, return longitude of 0 */ | |
11fdf7f2 TL |
171 | |
172 | /* calculate x,y coordinates with solution obtained */ | |
92f5a8d4 | 173 | if((D = acos(cos(lp_lat) * cos(C = 0.5 * lp_lon))) != 0.0) {/* Aitoff */ |
11fdf7f2 TL |
174 | x = 2. * D * cos(lp_lat) * sin(C) * (y = 1. / sin(D)); |
175 | y *= D * sin(lp_lat); | |
176 | } else | |
177 | x = y = 0.; | |
92f5a8d4 | 178 | if (this->m_proj_parm.mode == mode_winkel_tripel) { /* Winkel Tripel */ |
11fdf7f2 TL |
179 | x = (x + lp_lon * this->m_proj_parm.cosphi1) * 0.5; |
180 | y = (y + lp_lat) * 0.5; | |
181 | } | |
182 | /* if too far from given values of x,y, repeat with better approximation of phi,lam */ | |
92f5a8d4 | 183 | } while (((fabs(xy_x-x) > epsilon) || (fabs(xy_y-y) > epsilon)) && (round++ < max_round)); |
11fdf7f2 | 184 | |
92f5a8d4 TL |
185 | if (iter == max_iter && round == max_round) |
186 | { | |
187 | BOOST_THROW_EXCEPTION( projection_exception(error_non_convergent) ); | |
188 | //fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl); | |
189 | } | |
11fdf7f2 TL |
190 | } |
191 | ||
192 | static inline std::string get_name() | |
193 | { | |
194 | return "aitoff_spheroid"; | |
195 | } | |
196 | ||
197 | }; | |
198 | ||
92f5a8d4 TL |
199 | template <typename Parameters> |
200 | inline void setup(Parameters& par) | |
11fdf7f2 | 201 | { |
11fdf7f2 TL |
202 | par.es = 0.; |
203 | } | |
204 | ||
205 | ||
206 | // Aitoff | |
207 | template <typename Parameters, typename T> | |
208 | inline void setup_aitoff(Parameters& par, par_aitoff<T>& proj_parm) | |
209 | { | |
92f5a8d4 TL |
210 | proj_parm.mode = mode_aitoff; |
211 | setup(par); | |
11fdf7f2 TL |
212 | } |
213 | ||
214 | // Winkel Tripel | |
92f5a8d4 TL |
215 | template <typename Params, typename Parameters, typename T> |
216 | inline void setup_wintri(Params& params, Parameters& par, par_aitoff<T>& proj_parm) | |
11fdf7f2 | 217 | { |
92f5a8d4 TL |
218 | static const T two_div_pi = detail::two_div_pi<T>(); |
219 | ||
220 | T phi1; | |
11fdf7f2 | 221 | |
92f5a8d4 TL |
222 | proj_parm.mode = mode_winkel_tripel; |
223 | if (pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, phi1)) { | |
224 | if ((proj_parm.cosphi1 = cos(phi1)) == 0.) | |
225 | BOOST_THROW_EXCEPTION( projection_exception(error_lat_larger_than_90) ); | |
11fdf7f2 | 226 | } else /* 50d28' or phi1=acos(2/pi) */ |
92f5a8d4 TL |
227 | proj_parm.cosphi1 = two_div_pi; |
228 | setup(par); | |
11fdf7f2 TL |
229 | } |
230 | ||
231 | }} // namespace detail::aitoff | |
232 | #endif // doxygen | |
233 | ||
234 | /*! | |
235 | \brief Aitoff projection | |
236 | \ingroup projections | |
237 | \tparam Geographic latlong point type | |
238 | \tparam Cartesian xy point type | |
239 | \tparam Parameters parameter type | |
240 | \par Projection characteristics | |
241 | - Miscellaneous | |
242 | - Spheroid | |
243 | \par Example | |
244 | \image html ex_aitoff.gif | |
245 | */ | |
92f5a8d4 TL |
246 | template <typename T, typename Parameters> |
247 | struct aitoff_spheroid : public detail::aitoff::base_aitoff_spheroid<T, Parameters> | |
11fdf7f2 | 248 | { |
92f5a8d4 TL |
249 | template <typename Params> |
250 | inline aitoff_spheroid(Params const& , Parameters & par) | |
11fdf7f2 | 251 | { |
92f5a8d4 | 252 | detail::aitoff::setup_aitoff(par, this->m_proj_parm); |
11fdf7f2 TL |
253 | } |
254 | }; | |
255 | ||
256 | /*! | |
257 | \brief Winkel Tripel projection | |
258 | \ingroup projections | |
259 | \tparam Geographic latlong point type | |
260 | \tparam Cartesian xy point type | |
261 | \tparam Parameters parameter type | |
262 | \par Projection characteristics | |
263 | - Miscellaneous | |
264 | - Spheroid | |
265 | \par Projection parameters | |
266 | - lat_1: Latitude of first standard parallel (degrees) | |
267 | \par Example | |
268 | \image html ex_wintri.gif | |
269 | */ | |
92f5a8d4 TL |
270 | template <typename T, typename Parameters> |
271 | struct wintri_spheroid : public detail::aitoff::base_aitoff_spheroid<T, Parameters> | |
11fdf7f2 | 272 | { |
92f5a8d4 TL |
273 | template <typename Params> |
274 | inline wintri_spheroid(Params const& params, Parameters & par) | |
11fdf7f2 | 275 | { |
92f5a8d4 | 276 | detail::aitoff::setup_wintri(params, par, this->m_proj_parm); |
11fdf7f2 TL |
277 | } |
278 | }; | |
279 | ||
280 | #ifndef DOXYGEN_NO_DETAIL | |
281 | namespace detail | |
282 | { | |
283 | ||
284 | // Static projection | |
92f5a8d4 TL |
285 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_aitoff, aitoff_spheroid) |
286 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_wintri, wintri_spheroid) | |
11fdf7f2 TL |
287 | |
288 | // Factory entry(s) | |
92f5a8d4 TL |
289 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(aitoff_entry, aitoff_spheroid) |
290 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(wintri_entry, wintri_spheroid) | |
11fdf7f2 | 291 | |
92f5a8d4 | 292 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(aitoff_init) |
11fdf7f2 | 293 | { |
92f5a8d4 TL |
294 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(aitoff, aitoff_entry) |
295 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(wintri, wintri_entry) | |
11fdf7f2 TL |
296 | } |
297 | ||
298 | } // namespace detail | |
299 | #endif // doxygen | |
300 | ||
301 | } // namespace projections | |
302 | ||
303 | }} // namespace boost::geometry | |
304 | ||
305 | #endif // BOOST_GEOMETRY_PROJECTIONS_AITOFF_HPP | |
306 |