boost/geometry/srs/projections/proj/igh.hpp

   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, 2019.
   6 // Modifications copyright (c) 2017-2019, 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/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
  50 #include <boost/geometry/util/math.hpp>
  51
  52 #include <boost/shared_ptr.hpp>
  53
  54 namespace boost { namespace geometry
  55 {
  56
  57 namespace projections
  58 {
  59     #ifndef DOXYGEN_NO_DETAIL
  60     namespace detail { namespace igh
  61     {
  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.
  73
  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>
  82             struct par_igh
  83             {
  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                 }
 125             };
 126
 127             /* 40d 44' 11.8" [degrees] */
 128             template <typename T>
 129             inline T d4044118() { return (T(40) + T(44)/T(60.) + T(11.8)/T(3600.)) * geometry::math::d2r<T>(); }
 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
 156             static const double epsilon = 1.e-10; // allow a little 'slack' on zone edge positions
 157
 158             template <typename T, typename Parameters>
 159             struct base_igh_spheroid
 160             {
 161                 par_igh<T, Parameters> m_proj_parm;
 162
 163                 template <typename Params>
 164                 inline base_igh_spheroid(Params const& params, Parameters & par)
 165                     : m_proj_parm(params, par)
 166                 {}
 167
 168                 // FORWARD(s_forward)  spheroid
 169                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 170                 inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 171                 {
 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>();
 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
 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;
 202                 }
 203
 204                 // INVERSE(s_inverse)  spheroid
 205                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 206                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 207                 {
 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;
 221
 222                     const T y90 = this->m_proj_parm.dy0 + sqrt(c2); // lt=90 corresponds to y=y0+sqrt(2.0)
 223
 224                         int z = 0;
 225                         if (xy_y > y90+epsilon || xy_y < -y90+epsilon) // 0
 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
 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;
 252
 253                           switch (z) {
 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;
 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
 289             template <typename Params, typename Parameters, typename T>
 290             inline void setup_igh(Params const& params, Parameters& par, par_igh<T, Parameters>& proj_parm)
 291             {
 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>();
 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             */
 320
 321                     T lp_lam = 0, lp_phi = d4044118;
 322                     T xy1_x, xy1_y;
 323                     T xy3_x, xy3_y;
 324
 325                     // sinusoidal zones
 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);
 332
 333                     // mollweide zones
 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.
 339
 340                     // y0 ?
 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
 343                     // y0 + xy1_y = xy3_y for lt = 40d44'11.8"
 344                     proj_parm.dy0 = xy3_y - xy1_y;
 345
 346                     proj_parm.zones[0].y0 = proj_parm.dy0; // zone 1
 347
 348                     // mollweide zones (cont'd)
 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.;
 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     */
 374     template <typename T, typename Parameters>
 375     struct igh_spheroid : public detail::igh::base_igh_spheroid<T, Parameters>
 376     {
 377         template <typename Params>
 378         inline igh_spheroid(Params const& params, Parameters & par)
 379             : detail::igh::base_igh_spheroid<T, Parameters>(params, par)
 380         {
 381             detail::igh::setup_igh(params, par, this->m_proj_parm);
 382         }
 383     };
 384
 385     #ifndef DOXYGEN_NO_DETAIL
 386     namespace detail
 387     {
 388
 389         // Static projection
 390         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_igh, igh_spheroid)
 391
 392         // Factory entry(s)
 393         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(igh_entry, igh_spheroid)
 394
 395         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(igh_init)
 396         {
 397             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(igh, igh_entry)
 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