]>
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 | // Permission is hereby granted, free of charge, to any person obtaining a | |
23 | // copy of this software and associated documentation files (the "Software"), | |
24 | // to deal in the Software without restriction, including without limitation | |
25 | // the rights to use, copy, modify, merge, publish, distribute, sublicense, | |
26 | // and/or sell copies of the Software, and to permit persons to whom the | |
27 | // Software is furnished to do so, subject to the following conditions: | |
28 | ||
29 | // The above copyright notice and this permission notice shall be included | |
30 | // in all copies or substantial portions of the Software. | |
31 | ||
32 | // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS | |
33 | // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | |
34 | // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL | |
35 | // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | |
36 | // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING | |
37 | // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER | |
38 | // DEALINGS IN THE SOFTWARE. | |
39 | ||
92f5a8d4 TL |
40 | #ifndef BOOST_GEOMETRY_PROJECTIONS_IGH_HPP |
41 | #define BOOST_GEOMETRY_PROJECTIONS_IGH_HPP | |
11fdf7f2 TL |
42 | |
43 | #include <boost/geometry/srs/projections/impl/base_static.hpp> | |
44 | #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> | |
45 | #include <boost/geometry/srs/projections/impl/projects.hpp> | |
46 | #include <boost/geometry/srs/projections/impl/factory_entry.hpp> | |
47 | #include <boost/geometry/srs/projections/proj/gn_sinu.hpp> | |
48 | #include <boost/geometry/srs/projections/proj/moll.hpp> | |
49 | ||
92f5a8d4 | 50 | #include <boost/geometry/util/math.hpp> |
11fdf7f2 | 51 | |
92f5a8d4 | 52 | #include <boost/shared_ptr.hpp> |
11fdf7f2 | 53 | |
92f5a8d4 TL |
54 | namespace boost { namespace geometry |
55 | { | |
11fdf7f2 TL |
56 | |
57 | namespace projections | |
58 | { | |
59 | #ifndef DOXYGEN_NO_DETAIL | |
60 | namespace detail { namespace igh | |
61 | { | |
92f5a8d4 TL |
62 | template <typename T> |
63 | struct par_igh_zone | |
64 | { | |
65 | T x0; | |
66 | T y0; | |
67 | T lam0; | |
68 | }; | |
69 | ||
70 | // NOTE: x0, y0, lam0 are not used in moll nor sinu projections | |
71 | // so it is a waste of memory to keep 12 copies of projections | |
72 | // with parameters as in the original Proj4. | |
11fdf7f2 | 73 | |
92f5a8d4 TL |
74 | // TODO: It would be possible to further decrease the size of par_igh |
75 | // because spherical sinu and moll has constant parameters. | |
76 | // TODO: Furthermore there is no need to store par_igh_zone parameters | |
77 | // since they are constant for zones. In both fwd() and inv() there are | |
78 | // parts of code dependent on specific zones (if statements) anyway | |
79 | // so these parameters could be hardcoded there instead of stored. | |
80 | ||
81 | template <typename T, typename Parameters> | |
11fdf7f2 TL |
82 | struct par_igh |
83 | { | |
92f5a8d4 TL |
84 | moll_spheroid<T, Parameters> moll; |
85 | sinu_spheroid<T, Parameters> sinu; | |
86 | par_igh_zone<T> zones[12]; | |
87 | T dy0; | |
88 | ||
89 | // NOTE: The constructors of moll and sinu projections sets | |
90 | // par.es = 0 | |
91 | ||
92 | template <typename Params> | |
93 | inline par_igh(Params const& params, Parameters & par) | |
94 | : moll(params, par) | |
95 | , sinu(params, par) | |
96 | {} | |
97 | ||
98 | inline void fwd(int zone, Parameters const& par, T const& lp_lon, T const& lp_lat, T & xy_x, T & xy_y) const | |
99 | { | |
100 | if (zone <= 2 || zone >= 9) // 1, 2, 9, 10, 11, 12 | |
101 | moll.fwd(par, lp_lon, lp_lat, xy_x, xy_y); | |
102 | else // 3, 4, 5, 6, 7, 8 | |
103 | sinu.fwd(par, lp_lon, lp_lat, xy_x, xy_y); | |
104 | } | |
105 | ||
106 | inline void inv(int zone, Parameters const& par, T const& xy_x, T const& xy_y, T & lp_lon, T & lp_lat) const | |
107 | { | |
108 | if (zone <= 2 || zone >= 9) // 1, 2, 9, 10, 11, 12 | |
109 | moll.inv(par, xy_x, xy_y, lp_lon, lp_lat); | |
110 | else // 3, 4, 5, 6, 7, 8 | |
111 | sinu.inv(par, xy_x, xy_y, lp_lon, lp_lat); | |
112 | } | |
113 | ||
114 | inline void set_zone(int zone, T const& x_0, T const& y_0, T const& lon_0) | |
115 | { | |
116 | zones[zone - 1].x0 = x_0; | |
117 | zones[zone - 1].y0 = y_0; | |
118 | zones[zone - 1].lam0 = lon_0; | |
119 | } | |
120 | ||
121 | inline par_igh_zone<T> const& get_zone(int zone) const | |
122 | { | |
123 | return zones[zone - 1]; | |
124 | } | |
11fdf7f2 TL |
125 | }; |
126 | ||
92f5a8d4 | 127 | /* 40d 44' 11.8" [degrees] */ |
11fdf7f2 | 128 | template <typename T> |
92f5a8d4 | 129 | inline T d4044118() { return (T(40) + T(44)/T(60.) + T(11.8)/T(3600.)) * geometry::math::d2r<T>(); } |
11fdf7f2 TL |
130 | |
131 | template <typename T> | |
132 | inline T d10() { return T(10) * geometry::math::d2r<T>(); } | |
133 | template <typename T> | |
134 | inline T d20() { return T(20) * geometry::math::d2r<T>(); } | |
135 | template <typename T> | |
136 | inline T d30() { return T(30) * geometry::math::d2r<T>(); } | |
137 | template <typename T> | |
138 | inline T d40() { return T(40) * geometry::math::d2r<T>(); } | |
139 | template <typename T> | |
140 | inline T d50() { return T(50) * geometry::math::d2r<T>(); } | |
141 | template <typename T> | |
142 | inline T d60() { return T(60) * geometry::math::d2r<T>(); } | |
143 | template <typename T> | |
144 | inline T d80() { return T(80) * geometry::math::d2r<T>(); } | |
145 | template <typename T> | |
146 | inline T d90() { return T(90) * geometry::math::d2r<T>(); } | |
147 | template <typename T> | |
148 | inline T d100() { return T(100) * geometry::math::d2r<T>(); } | |
149 | template <typename T> | |
150 | inline T d140() { return T(140) * geometry::math::d2r<T>(); } | |
151 | template <typename T> | |
152 | inline T d160() { return T(160) * geometry::math::d2r<T>(); } | |
153 | template <typename T> | |
154 | inline T d180() { return T(180) * geometry::math::d2r<T>(); } | |
155 | ||
92f5a8d4 | 156 | static const double epsilon = 1.e-10; // allow a little 'slack' on zone edge positions |
11fdf7f2 | 157 | |
92f5a8d4 TL |
158 | template <typename T, typename Parameters> |
159 | struct base_igh_spheroid | |
11fdf7f2 | 160 | { |
92f5a8d4 | 161 | par_igh<T, Parameters> m_proj_parm; |
11fdf7f2 | 162 | |
92f5a8d4 TL |
163 | template <typename Params> |
164 | inline base_igh_spheroid(Params const& params, Parameters & par) | |
165 | : m_proj_parm(params, par) | |
166 | {} | |
11fdf7f2 TL |
167 | |
168 | // FORWARD(s_forward) spheroid | |
169 | // Project coordinates from geographic (lon, lat) to cartesian (x, y) | |
92f5a8d4 | 170 | inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const |
11fdf7f2 | 171 | { |
92f5a8d4 TL |
172 | static const T d4044118 = igh::d4044118<T>(); |
173 | static const T d20 = igh::d20<T>(); | |
174 | static const T d40 = igh::d40<T>(); | |
175 | static const T d80 = igh::d80<T>(); | |
176 | static const T d100 = igh::d100<T>(); | |
11fdf7f2 TL |
177 | |
178 | int z; | |
179 | if (lp_lat >= d4044118) { // 1|2 | |
180 | z = (lp_lon <= -d40 ? 1: 2); | |
181 | } | |
182 | else if (lp_lat >= 0) { // 3|4 | |
183 | z = (lp_lon <= -d40 ? 3: 4); | |
184 | } | |
185 | else if (lp_lat >= -d4044118) { // 5|6|7|8 | |
186 | if (lp_lon <= -d100) z = 5; // 5 | |
187 | else if (lp_lon <= -d20) z = 6; // 6 | |
188 | else if (lp_lon <= d80) z = 7; // 7 | |
189 | else z = 8; // 8 | |
190 | } | |
191 | else { // 9|10|11|12 | |
192 | if (lp_lon <= -d100) z = 9; // 9 | |
193 | else if (lp_lon <= -d20) z = 10; // 10 | |
194 | else if (lp_lon <= d80) z = 11; // 11 | |
195 | else z = 12; // 12 | |
196 | } | |
197 | ||
92f5a8d4 TL |
198 | lp_lon -= this->m_proj_parm.get_zone(z).lam0; |
199 | this->m_proj_parm.fwd(z, par, lp_lon, lp_lat, xy_x, xy_y); | |
200 | xy_x += this->m_proj_parm.get_zone(z).x0; | |
201 | xy_y += this->m_proj_parm.get_zone(z).y0; | |
11fdf7f2 TL |
202 | } |
203 | ||
204 | // INVERSE(s_inverse) spheroid | |
205 | // Project coordinates from cartesian (x, y) to geographic (lon, lat) | |
92f5a8d4 | 206 | inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const |
11fdf7f2 | 207 | { |
92f5a8d4 TL |
208 | static const T d4044118 = igh::d4044118<T>(); |
209 | static const T d10 = igh::d10<T>(); | |
210 | static const T d20 = igh::d20<T>(); | |
211 | static const T d40 = igh::d40<T>(); | |
212 | static const T d50 = igh::d50<T>(); | |
213 | static const T d60 = igh::d60<T>(); | |
214 | static const T d80 = igh::d80<T>(); | |
215 | static const T d90 = igh::d90<T>(); | |
216 | static const T d100 = igh::d100<T>(); | |
217 | static const T d160 = igh::d160<T>(); | |
218 | static const T d180 = igh::d180<T>(); | |
219 | ||
220 | static const T c2 = 2.0; | |
11fdf7f2 | 221 | |
92f5a8d4 | 222 | const T y90 = this->m_proj_parm.dy0 + sqrt(c2); // lt=90 corresponds to y=y0+sqrt(2.0) |
11fdf7f2 TL |
223 | |
224 | int z = 0; | |
92f5a8d4 | 225 | if (xy_y > y90+epsilon || xy_y < -y90+epsilon) // 0 |
11fdf7f2 TL |
226 | z = 0; |
227 | else if (xy_y >= d4044118) // 1|2 | |
228 | z = (xy_x <= -d40? 1: 2); | |
229 | else if (xy_y >= 0) // 3|4 | |
230 | z = (xy_x <= -d40? 3: 4); | |
231 | else if (xy_y >= -d4044118) { // 5|6|7|8 | |
232 | if (xy_x <= -d100) z = 5; // 5 | |
233 | else if (xy_x <= -d20) z = 6; // 6 | |
234 | else if (xy_x <= d80) z = 7; // 7 | |
235 | else z = 8; // 8 | |
236 | } | |
237 | else { // 9|10|11|12 | |
238 | if (xy_x <= -d100) z = 9; // 9 | |
239 | else if (xy_x <= -d20) z = 10; // 10 | |
240 | else if (xy_x <= d80) z = 11; // 11 | |
241 | else z = 12; // 12 | |
242 | } | |
243 | ||
244 | if (z) | |
245 | { | |
246 | int ok = 0; | |
247 | ||
92f5a8d4 TL |
248 | xy_x -= this->m_proj_parm.get_zone(z).x0; |
249 | xy_y -= this->m_proj_parm.get_zone(z).y0; | |
250 | this->m_proj_parm.inv(z, par, xy_x, xy_y, lp_lon, lp_lat); | |
251 | lp_lon += this->m_proj_parm.get_zone(z).lam0; | |
11fdf7f2 TL |
252 | |
253 | switch (z) { | |
92f5a8d4 TL |
254 | case 1: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d40+epsilon) || |
255 | ((lp_lon >= -d40-epsilon && lp_lon <= -d10+epsilon) && | |
256 | (lp_lat >= d60-epsilon && lp_lat <= d90+epsilon)); break; | |
257 | case 2: ok = (lp_lon >= -d40-epsilon && lp_lon <= d180+epsilon) || | |
258 | ((lp_lon >= -d180-epsilon && lp_lon <= -d160+epsilon) && | |
259 | (lp_lat >= d50-epsilon && lp_lat <= d90+epsilon)) || | |
260 | ((lp_lon >= -d50-epsilon && lp_lon <= -d40+epsilon) && | |
261 | (lp_lat >= d60-epsilon && lp_lat <= d90+epsilon)); break; | |
262 | case 3: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d40+epsilon); break; | |
263 | case 4: ok = (lp_lon >= -d40-epsilon && lp_lon <= d180+epsilon); break; | |
264 | case 5: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d100+epsilon); break; | |
265 | case 6: ok = (lp_lon >= -d100-epsilon && lp_lon <= -d20+epsilon); break; | |
266 | case 7: ok = (lp_lon >= -d20-epsilon && lp_lon <= d80+epsilon); break; | |
267 | case 8: ok = (lp_lon >= d80-epsilon && lp_lon <= d180+epsilon); break; | |
268 | case 9: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d100+epsilon); break; | |
269 | case 10: ok = (lp_lon >= -d100-epsilon && lp_lon <= -d20+epsilon); break; | |
270 | case 11: ok = (lp_lon >= -d20-epsilon && lp_lon <= d80+epsilon); break; | |
271 | case 12: ok = (lp_lon >= d80-epsilon && lp_lon <= d180+epsilon); break; | |
11fdf7f2 TL |
272 | } |
273 | ||
274 | z = (!ok? 0: z); // projectable? | |
275 | } | |
276 | // if (!z) pj_errno = -15; // invalid x or y | |
277 | if (!z) lp_lon = HUGE_VAL; | |
278 | if (!z) lp_lat = HUGE_VAL; | |
279 | } | |
280 | ||
281 | static inline std::string get_name() | |
282 | { | |
283 | return "igh_spheroid"; | |
284 | } | |
285 | ||
286 | }; | |
287 | ||
288 | // Interrupted Goode Homolosine | |
92f5a8d4 TL |
289 | template <typename Params, typename Parameters, typename T> |
290 | inline void setup_igh(Params const& params, Parameters& par, par_igh<T, Parameters>& proj_parm) | |
11fdf7f2 | 291 | { |
92f5a8d4 TL |
292 | static const T d0 = 0; |
293 | static const T d4044118 = igh::d4044118<T>(); | |
294 | static const T d20 = igh::d20<T>(); | |
295 | static const T d30 = igh::d30<T>(); | |
296 | static const T d60 = igh::d60<T>(); | |
297 | static const T d100 = igh::d100<T>(); | |
298 | static const T d140 = igh::d140<T>(); | |
299 | static const T d160 = igh::d160<T>(); | |
11fdf7f2 TL |
300 | |
301 | /* | |
302 | Zones: | |
303 | ||
304 | -180 -40 180 | |
305 | +--------------+-------------------------+ Zones 1,2,9,10,11 & 12: | |
306 | |1 |2 | Mollweide projection | |
307 | | | | | |
308 | +--------------+-------------------------+ Zones 3,4,5,6,7 & 8: | |
309 | |3 |4 | Sinusoidal projection | |
310 | | | | | |
311 | 0 +-------+------+-+-----------+-----------+ | |
312 | |5 |6 |7 |8 | | |
313 | | | | | | | |
314 | +-------+--------+-----------+-----------+ | |
315 | |9 |10 |11 |12 | | |
316 | | | | | | | |
317 | +-------+--------+-----------+-----------+ | |
318 | -180 -100 -20 80 180 | |
319 | */ | |
92f5a8d4 TL |
320 | |
321 | T lp_lam = 0, lp_phi = d4044118; | |
322 | T xy1_x, xy1_y; | |
323 | T xy3_x, xy3_y; | |
11fdf7f2 TL |
324 | |
325 | // sinusoidal zones | |
92f5a8d4 TL |
326 | proj_parm.set_zone(3, -d100, d0, -d100); |
327 | proj_parm.set_zone(4, d30, d0, d30); | |
328 | proj_parm.set_zone(5, -d160, d0, -d160); | |
329 | proj_parm.set_zone(6, -d60, d0, -d60); | |
330 | proj_parm.set_zone(7, d20, d0, d20); | |
331 | proj_parm.set_zone(8, d140, d0, d140); | |
11fdf7f2 TL |
332 | |
333 | // mollweide zones | |
92f5a8d4 TL |
334 | proj_parm.set_zone(1, -d100, d0, -d100); |
335 | ||
336 | // NOTE: x0, y0, lam0 are not used in moll nor sinu fwd | |
337 | // so the order of initialization doesn't matter that much. | |
338 | // But keep the original one from Proj4. | |
11fdf7f2 TL |
339 | |
340 | // y0 ? | |
92f5a8d4 TL |
341 | proj_parm.fwd(1, par, lp_lam, lp_phi, xy1_x, xy1_y); // zone 1 |
342 | proj_parm.fwd(3, par, lp_lam, lp_phi, xy3_x, xy3_y); // zone 3 | |
11fdf7f2 TL |
343 | // y0 + xy1_y = xy3_y for lt = 40d44'11.8" |
344 | proj_parm.dy0 = xy3_y - xy1_y; | |
345 | ||
92f5a8d4 | 346 | proj_parm.zones[0].y0 = proj_parm.dy0; // zone 1 |
11fdf7f2 TL |
347 | |
348 | // mollweide zones (cont'd) | |
92f5a8d4 TL |
349 | proj_parm.set_zone(2, d30, proj_parm.dy0, d30); |
350 | proj_parm.set_zone(9, -d160, -proj_parm.dy0, -d160); | |
351 | proj_parm.set_zone(10, -d60, -proj_parm.dy0, -d60); | |
352 | proj_parm.set_zone(11, d20, -proj_parm.dy0, d20); | |
353 | proj_parm.set_zone(12, d140, -proj_parm.dy0, d140); | |
354 | ||
355 | // NOTE: Already done before in sinu and moll constructor | |
356 | //par.es = 0.; | |
11fdf7f2 TL |
357 | } |
358 | ||
359 | }} // namespace detail::igh | |
360 | #endif // doxygen | |
361 | ||
362 | /*! | |
363 | \brief Interrupted Goode Homolosine projection | |
364 | \ingroup projections | |
365 | \tparam Geographic latlong point type | |
366 | \tparam Cartesian xy point type | |
367 | \tparam Parameters parameter type | |
368 | \par Projection characteristics | |
369 | - Pseudocylindrical | |
370 | - Spheroid | |
371 | \par Example | |
372 | \image html ex_igh.gif | |
373 | */ | |
92f5a8d4 TL |
374 | template <typename T, typename Parameters> |
375 | struct igh_spheroid : public detail::igh::base_igh_spheroid<T, Parameters> | |
11fdf7f2 | 376 | { |
92f5a8d4 TL |
377 | template <typename Params> |
378 | inline igh_spheroid(Params const& params, Parameters & par) | |
379 | : detail::igh::base_igh_spheroid<T, Parameters>(params, par) | |
11fdf7f2 | 380 | { |
92f5a8d4 | 381 | detail::igh::setup_igh(params, par, this->m_proj_parm); |
11fdf7f2 TL |
382 | } |
383 | }; | |
384 | ||
385 | #ifndef DOXYGEN_NO_DETAIL | |
386 | namespace detail | |
387 | { | |
388 | ||
389 | // Static projection | |
92f5a8d4 | 390 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_igh, igh_spheroid) |
11fdf7f2 TL |
391 | |
392 | // Factory entry(s) | |
92f5a8d4 | 393 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(igh_entry, igh_spheroid) |
11fdf7f2 | 394 | |
92f5a8d4 | 395 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(igh_init) |
11fdf7f2 | 396 | { |
92f5a8d4 | 397 | BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(igh, igh_entry) |
11fdf7f2 TL |
398 | } |
399 | ||
400 | } // namespace detail | |
401 | #endif // doxygen | |
402 | ||
403 | } // namespace projections | |
404 | ||
405 | }} // namespace boost::geometry | |
406 | ||
407 | #endif // BOOST_GEOMETRY_PROJECTIONS_IGH_HPP | |
408 |