1 // Boost.Geometry - gis-projections (based on PROJ4) 2 3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. 4 5 // This file was modified by Oracle on 2017, 2018. 6 // Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. 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 18 // Last updated version of proj: 5.0.0 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 40 #ifndef BOOST_GEOMETRY_PROJECTIONS_IGH_HPP 41 #define BOOST_GEOMETRY_PROJECTIONS_IGH_HPP 42 43 #include <boost/geometry/util/math.hpp> 44 #include <boost/shared_ptr.hpp> 45 46 #include <boost/geometry/srs/projections/impl/base_static.hpp> 47 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> 48 #include <boost/geometry/srs/projections/impl/projects.hpp> 49 #include <boost/geometry/srs/projections/impl/factory_entry.hpp> 50 #include <boost/geometry/srs/projections/proj/gn_sinu.hpp> 51 #include <boost/geometry/srs/projections/proj/moll.hpp> 52 53 namespace boost { namespace geometry 54 { 55 56 namespace projections 57 { 58 #ifndef DOXYGEN_NO_DETAIL 59 namespace detail { namespace igh 60 { 61 // TODO: consider replacing dynamically created projections 62 // with member objects 63 template <typename T, typename Parameters> 64 struct par_igh 65 { 66 boost::shared_ptr<base_v<T, Parameters> > pj[12]; 67 T dy0; 68 }; 69 70 /* 40d 44' 11.8" [degrees] */ 71 template <typename T> d4044118()72 inline T d4044118() { return (T(40) + T(44)/T(60.) + T(11.8)/T(3600.)) * geometry::math::d2r<T>(); } 73 74 template <typename T> d10()75 inline T d10() { return T(10) * geometry::math::d2r<T>(); } 76 template <typename T> d20()77 inline T d20() { return T(20) * geometry::math::d2r<T>(); } 78 template <typename T> d30()79 inline T d30() { return T(30) * geometry::math::d2r<T>(); } 80 template <typename T> d40()81 inline T d40() { return T(40) * geometry::math::d2r<T>(); } 82 template <typename T> d50()83 inline T d50() { return T(50) * geometry::math::d2r<T>(); } 84 template <typename T> d60()85 inline T d60() { return T(60) * geometry::math::d2r<T>(); } 86 template <typename T> d80()87 inline T d80() { return T(80) * geometry::math::d2r<T>(); } 88 template <typename T> d90()89 inline T d90() { return T(90) * geometry::math::d2r<T>(); } 90 template <typename T> d100()91 inline T d100() { return T(100) * geometry::math::d2r<T>(); } 92 template <typename T> d140()93 inline T d140() { return T(140) * geometry::math::d2r<T>(); } 94 template <typename T> d160()95 inline T d160() { return T(160) * geometry::math::d2r<T>(); } 96 template <typename T> d180()97 inline T d180() { return T(180) * geometry::math::d2r<T>(); } 98 99 static const double epsilon = 1.e-10; // allow a little 'slack' on zone edge positions 100 101 // Converted from #define SETUP(n, proj, x_0, y_0, lon_0) 102 template <template <typename, typename, typename> class Entry, typename Params, typename Parameters, typename T> do_setup(int n,Params const & params,Parameters const & par,par_igh<T,Parameters> & proj_parm,T const & x_0,T const & y_0,T const & lon_0)103 inline void do_setup(int n, Params const& params, Parameters const& par, par_igh<T, Parameters>& proj_parm, 104 T const& x_0, T const& y_0, 105 T const& lon_0) 106 { 107 // NOTE: in the original proj4 these projections are initialized 108 // with zeroed parameters which could be done here as well instead 109 // of initializing with parent projection's parameters. 110 Entry<Params, T, Parameters> entry; 111 proj_parm.pj[n-1].reset(entry.create_new(params, par)); 112 proj_parm.pj[n-1]->mutable_params().x0 = x_0; 113 proj_parm.pj[n-1]->mutable_params().y0 = y_0; 114 proj_parm.pj[n-1]->mutable_params().lam0 = lon_0; 115 } 116 117 // template class, using CRTP to implement forward/inverse 118 template <typename T, typename Parameters> 119 struct base_igh_spheroid 120 : public base_t_fi<base_igh_spheroid<T, Parameters>, T, Parameters> 121 { 122 par_igh<T, Parameters> m_proj_parm; 123 base_igh_spheroidboost::geometry::projections::detail::igh::base_igh_spheroid124 inline base_igh_spheroid(const Parameters& par) 125 : base_t_fi<base_igh_spheroid<T, Parameters>, T, Parameters>(*this, par) 126 {} 127 128 // FORWARD(s_forward) spheroid 129 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::igh::base_igh_spheroid130 inline void fwd(T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 131 { 132 static const T d4044118 = igh::d4044118<T>(); 133 static const T d20 = igh::d20<T>(); 134 static const T d40 = igh::d40<T>(); 135 static const T d80 = igh::d80<T>(); 136 static const T d100 = igh::d100<T>(); 137 138 int z; 139 if (lp_lat >= d4044118) { // 1|2 140 z = (lp_lon <= -d40 ? 1: 2); 141 } 142 else if (lp_lat >= 0) { // 3|4 143 z = (lp_lon <= -d40 ? 3: 4); 144 } 145 else if (lp_lat >= -d4044118) { // 5|6|7|8 146 if (lp_lon <= -d100) z = 5; // 5 147 else if (lp_lon <= -d20) z = 6; // 6 148 else if (lp_lon <= d80) z = 7; // 7 149 else z = 8; // 8 150 } 151 else { // 9|10|11|12 152 if (lp_lon <= -d100) z = 9; // 9 153 else if (lp_lon <= -d20) z = 10; // 10 154 else if (lp_lon <= d80) z = 11; // 11 155 else z = 12; // 12 156 } 157 158 lp_lon -= this->m_proj_parm.pj[z-1]->params().lam0; 159 this->m_proj_parm.pj[z-1]->fwd(lp_lon, lp_lat, xy_x, xy_y); 160 xy_x += this->m_proj_parm.pj[z-1]->params().x0; 161 xy_y += this->m_proj_parm.pj[z-1]->params().y0; 162 } 163 164 // INVERSE(s_inverse) spheroid 165 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::igh::base_igh_spheroid166 inline void inv(T xy_x, T xy_y, T& lp_lon, T& lp_lat) const 167 { 168 static const T d4044118 = igh::d4044118<T>(); 169 static const T d10 = igh::d10<T>(); 170 static const T d20 = igh::d20<T>(); 171 static const T d40 = igh::d40<T>(); 172 static const T d50 = igh::d50<T>(); 173 static const T d60 = igh::d60<T>(); 174 static const T d80 = igh::d80<T>(); 175 static const T d90 = igh::d90<T>(); 176 static const T d100 = igh::d100<T>(); 177 static const T d160 = igh::d160<T>(); 178 static const T d180 = igh::d180<T>(); 179 180 static const T c2 = 2.0; 181 182 const T y90 = this->m_proj_parm.dy0 + sqrt(c2); // lt=90 corresponds to y=y0+sqrt(2.0) 183 184 int z = 0; 185 if (xy_y > y90+epsilon || xy_y < -y90+epsilon) // 0 186 z = 0; 187 else if (xy_y >= d4044118) // 1|2 188 z = (xy_x <= -d40? 1: 2); 189 else if (xy_y >= 0) // 3|4 190 z = (xy_x <= -d40? 3: 4); 191 else if (xy_y >= -d4044118) { // 5|6|7|8 192 if (xy_x <= -d100) z = 5; // 5 193 else if (xy_x <= -d20) z = 6; // 6 194 else if (xy_x <= d80) z = 7; // 7 195 else z = 8; // 8 196 } 197 else { // 9|10|11|12 198 if (xy_x <= -d100) z = 9; // 9 199 else if (xy_x <= -d20) z = 10; // 10 200 else if (xy_x <= d80) z = 11; // 11 201 else z = 12; // 12 202 } 203 204 if (z) 205 { 206 int ok = 0; 207 208 xy_x -= this->m_proj_parm.pj[z-1]->params().x0; 209 xy_y -= this->m_proj_parm.pj[z-1]->params().y0; 210 this->m_proj_parm.pj[z-1]->inv(xy_x, xy_y, lp_lon, lp_lat); 211 lp_lon += this->m_proj_parm.pj[z-1]->params().lam0; 212 213 switch (z) { 214 case 1: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d40+epsilon) || 215 ((lp_lon >= -d40-epsilon && lp_lon <= -d10+epsilon) && 216 (lp_lat >= d60-epsilon && lp_lat <= d90+epsilon)); break; 217 case 2: ok = (lp_lon >= -d40-epsilon && lp_lon <= d180+epsilon) || 218 ((lp_lon >= -d180-epsilon && lp_lon <= -d160+epsilon) && 219 (lp_lat >= d50-epsilon && lp_lat <= d90+epsilon)) || 220 ((lp_lon >= -d50-epsilon && lp_lon <= -d40+epsilon) && 221 (lp_lat >= d60-epsilon && lp_lat <= d90+epsilon)); break; 222 case 3: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d40+epsilon); break; 223 case 4: ok = (lp_lon >= -d40-epsilon && lp_lon <= d180+epsilon); break; 224 case 5: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d100+epsilon); break; 225 case 6: ok = (lp_lon >= -d100-epsilon && lp_lon <= -d20+epsilon); break; 226 case 7: ok = (lp_lon >= -d20-epsilon && lp_lon <= d80+epsilon); break; 227 case 8: ok = (lp_lon >= d80-epsilon && lp_lon <= d180+epsilon); break; 228 case 9: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d100+epsilon); break; 229 case 10: ok = (lp_lon >= -d100-epsilon && lp_lon <= -d20+epsilon); break; 230 case 11: ok = (lp_lon >= -d20-epsilon && lp_lon <= d80+epsilon); break; 231 case 12: ok = (lp_lon >= d80-epsilon && lp_lon <= d180+epsilon); break; 232 } 233 234 z = (!ok? 0: z); // projectable? 235 } 236 // if (!z) pj_errno = -15; // invalid x or y 237 if (!z) lp_lon = HUGE_VAL; 238 if (!z) lp_lat = HUGE_VAL; 239 } 240 get_nameboost::geometry::projections::detail::igh::base_igh_spheroid241 static inline std::string get_name() 242 { 243 return "igh_spheroid"; 244 } 245 246 }; 247 248 // Interrupted Goode Homolosine 249 template <typename Params, typename Parameters, typename T> setup_igh(Params const & params,Parameters & par,par_igh<T,Parameters> & proj_parm)250 inline void setup_igh(Params const& params, Parameters& par, par_igh<T, Parameters>& proj_parm) 251 { 252 static const T d0 = 0; 253 static const T d4044118 = igh::d4044118<T>(); 254 static const T d20 = igh::d20<T>(); 255 static const T d30 = igh::d30<T>(); 256 static const T d60 = igh::d60<T>(); 257 static const T d100 = igh::d100<T>(); 258 static const T d140 = igh::d140<T>(); 259 static const T d160 = igh::d160<T>(); 260 261 /* 262 Zones: 263 264 -180 -40 180 265 +--------------+-------------------------+ Zones 1,2,9,10,11 & 12: 266 |1 |2 | Mollweide projection 267 | | | 268 +--------------+-------------------------+ Zones 3,4,5,6,7 & 8: 269 |3 |4 | Sinusoidal projection 270 | | | 271 0 +-------+------+-+-----------+-----------+ 272 |5 |6 |7 |8 | 273 | | | | | 274 +-------+--------+-----------+-----------+ 275 |9 |10 |11 |12 | 276 | | | | | 277 +-------+--------+-----------+-----------+ 278 -180 -100 -20 80 180 279 */ 280 281 T lp_lam = 0, lp_phi = d4044118; 282 T xy1_x, xy1_y; 283 T xy3_x, xy3_y; 284 285 // IMPORTANT: Force spherical sinu projection 286 // This is required because unlike in the original proj4 here 287 // parameters are used to initialize underlying projections. 288 // In the original code zeroed parameters are passed which 289 // could be done here as well though. 290 par.es = 0.; 291 292 // sinusoidal zones 293 do_setup<sinu_entry>(3, params, par, proj_parm, -d100, d0, -d100); 294 do_setup<sinu_entry>(4, params, par, proj_parm, d30, d0, d30); 295 do_setup<sinu_entry>(5, params, par, proj_parm, -d160, d0, -d160); 296 do_setup<sinu_entry>(6, params, par, proj_parm, -d60, d0, -d60); 297 do_setup<sinu_entry>(7, params, par, proj_parm, d20, d0, d20); 298 do_setup<sinu_entry>(8, params, par, proj_parm, d140, d0, d140); 299 300 // mollweide zones 301 do_setup<moll_entry>(1, params, par, proj_parm, -d100, d0, -d100); 302 303 // y0 ? 304 proj_parm.pj[0]->fwd(lp_lam, lp_phi, xy1_x, xy1_y); // zone 1 305 proj_parm.pj[2]->fwd(lp_lam, lp_phi, xy3_x, xy3_y); // zone 3 306 // y0 + xy1_y = xy3_y for lt = 40d44'11.8" 307 proj_parm.dy0 = xy3_y - xy1_y; 308 309 proj_parm.pj[0]->mutable_params().y0 = proj_parm.dy0; 310 311 // mollweide zones (cont'd) 312 do_setup<moll_entry>( 2, params, par, proj_parm, d30, proj_parm.dy0, d30); 313 do_setup<moll_entry>( 9, params, par, proj_parm, -d160, -proj_parm.dy0, -d160); 314 do_setup<moll_entry>(10, params, par, proj_parm, -d60, -proj_parm.dy0, -d60); 315 do_setup<moll_entry>(11, params, par, proj_parm, d20, -proj_parm.dy0, d20); 316 do_setup<moll_entry>(12, params, par, proj_parm, d140, -proj_parm.dy0, d140); 317 318 // Already done before 319 //par.es = 0.; 320 } 321 322 }} // namespace detail::igh 323 #endif // doxygen 324 325 /*! 326 \brief Interrupted Goode Homolosine projection 327 \ingroup projections 328 \tparam Geographic latlong point type 329 \tparam Cartesian xy point type 330 \tparam Parameters parameter type 331 \par Projection characteristics 332 - Pseudocylindrical 333 - Spheroid 334 \par Example 335 \image html ex_igh.gif 336 */ 337 template <typename T, typename Parameters> 338 struct igh_spheroid : public detail::igh::base_igh_spheroid<T, Parameters> 339 { 340 template <typename Params> igh_spheroidboost::geometry::projections::igh_spheroid341 inline igh_spheroid(Params const& params, Parameters const& par) 342 : detail::igh::base_igh_spheroid<T, Parameters>(par) 343 { 344 detail::igh::setup_igh(params, this->m_par, this->m_proj_parm); 345 } 346 }; 347 348 #ifndef DOXYGEN_NO_DETAIL 349 namespace detail 350 { 351 352 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::spar::proj_igh,igh_spheroid,igh_spheroid)353 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::spar::proj_igh, igh_spheroid, igh_spheroid) 354 355 // Factory entry(s) 356 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(igh_entry, igh_spheroid) 357 358 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(igh_init) 359 { 360 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(igh, igh_entry) 361 } 362 363 } // namespace detail 364 #endif // doxygen 365 366 } // namespace projections 367 368 }} // namespace boost::geometry 369 370 #endif // BOOST_GEOMETRY_PROJECTIONS_IGH_HPP 371 372