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

   1 #ifndef BOOST_GEOMETRY_PROJECTIONS_IGH_HPP
   2 #define BOOST_GEOMETRY_PROJECTIONS_IGH_HPP
   3
   4 // Boost.Geometry - extensions-gis-projections (based on PROJ4)
   5 // This file is automatically generated. DO NOT EDIT.
   6
   7 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
   8
   9 // This file was modified by Oracle on 2017.
  10 // Modifications copyright (c) 2017, Oracle and/or its affiliates.
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
  12
  13 // Use, modification and distribution is subject to the Boost Software License,
  14 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  15 // http://www.boost.org/LICENSE_1_0.txt)
  16
  17 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  18 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  19 // PROJ4 is maintained by Frank Warmerdam
  20 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
  21
  22 // Last updated version of proj: 4.9.1
  23
  24 // Original copyright notice:
  25
  26 // Permission is hereby granted, free of charge, to any person obtaining a
  27 // copy of this software and associated documentation files (the "Software"),
  28 // to deal in the Software without restriction, including without limitation
  29 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  30 // and/or sell copies of the Software, and to permit persons to whom the
  31 // Software is furnished to do so, subject to the following conditions:
  32
  33 // The above copyright notice and this permission notice shall be included
  34 // in all copies or substantial portions of the Software.
  35
  36 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  37 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  38 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  39 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  40 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  41 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  42 // DEALINGS IN THE SOFTWARE.
  43
  44 #include <boost/geometry/util/math.hpp>
  45 #include <boost/shared_ptr.hpp>
  46
  47 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  48 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  49 #include <boost/geometry/srs/projections/impl/projects.hpp>
  50 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  51 #include <boost/geometry/srs/projections/proj/gn_sinu.hpp>
  52 #include <boost/geometry/srs/projections/proj/moll.hpp>
  53
  54 namespace boost { namespace geometry
  55 {
  56
  57 namespace srs { namespace par4
  58 {
  59     struct igh {};
  60
  61 }} //namespace srs::par4
  62
  63 namespace projections
  64 {
  65     #ifndef DOXYGEN_NO_DETAIL
  66     namespace detail { namespace igh
  67     {
  68
  69             template <typename CalculationType, typename Parameters>
  70             struct par_igh
  71             {
  72                 boost::shared_ptr<base_v<CalculationType, Parameters> > pj[12];
  73                 CalculationType dy0;
  74             };
  75
  76             template <typename T>
  77             inline T d4044118() { return (T(40) + T(44)/T(60.) + T(11.8)/T(3600.)) * geometry::math::d2r<T>(); } // 40d 44' 11.8" [degrees]
  78
  79             template <typename T>
  80             inline T d10() { return T(10) * geometry::math::d2r<T>(); }
  81             template <typename T>
  82             inline T d20() { return T(20) * geometry::math::d2r<T>(); }
  83             template <typename T>
  84             inline T d30() { return T(30) * geometry::math::d2r<T>(); }
  85             template <typename T>
  86             inline T d40() { return T(40) * geometry::math::d2r<T>(); }
  87             template <typename T>
  88             inline T d50() { return T(50) * geometry::math::d2r<T>(); }
  89             template <typename T>
  90             inline T d60() { return T(60) * geometry::math::d2r<T>(); }
  91             template <typename T>
  92             inline T d80() { return T(80) * geometry::math::d2r<T>(); }
  93             template <typename T>
  94             inline T d90() { return T(90) * geometry::math::d2r<T>(); }
  95             template <typename T>
  96             inline T d100() { return T(100) * geometry::math::d2r<T>(); }
  97             template <typename T>
  98             inline T d140() { return T(140) * geometry::math::d2r<T>(); }
  99             template <typename T>
 100             inline T d160() { return T(160) * geometry::math::d2r<T>(); }
 101             template <typename T>
 102             inline T d180() { return T(180) * geometry::math::d2r<T>(); }
 103
 104             static const double EPSLN = 1.e-10; // allow a little 'slack' on zone edge positions
 105
 106             // Converted from #define SETUP(n, proj, x_0, y_0, lon_0)
 107             template <template <typename, typename> class Entry, typename Parameters, typename CalculationType>
 108             inline void do_setup(int n, Parameters const& par, par_igh<CalculationType, Parameters>& proj_parm,
 109                                  CalculationType const& x_0, CalculationType const& y_0,
 110                                  CalculationType const& lon_0)
 111             {
 112                 Entry<CalculationType, Parameters> entry;
 113                 proj_parm.pj[n-1].reset(entry.create_new(par));
 114                 proj_parm.pj[n-1]->mutable_params().x0 = x_0;
 115                 proj_parm.pj[n-1]->mutable_params().y0 = y_0;
 116                 proj_parm.pj[n-1]->mutable_params().lam0 = lon_0;
 117             }
 118
 119             // template class, using CRTP to implement forward/inverse
 120             template <typename CalculationType, typename Parameters>
 121             struct base_igh_spheroid : public base_t_fi<base_igh_spheroid<CalculationType, Parameters>,
 122                      CalculationType, Parameters>
 123             {
 124
 125                 typedef CalculationType geographic_type;
 126                 typedef CalculationType cartesian_type;
 127
 128                 par_igh<CalculationType, Parameters> m_proj_parm;
 129
 130                 inline base_igh_spheroid(const Parameters& par)
 131                     : base_t_fi<base_igh_spheroid<CalculationType, Parameters>,
 132                      CalculationType, Parameters>(*this, par) {}
 133
 134                 // FORWARD(s_forward)  spheroid
 135                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 136                 inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
 137                 {
 138                     static const CalculationType d4044118 = igh::d4044118<CalculationType>();
 139                     static const CalculationType d20  =  igh::d20<CalculationType>();
 140                     static const CalculationType d40  =  igh::d40<CalculationType>();
 141                     static const CalculationType d80  =  igh::d80<CalculationType>();
 142                     static const CalculationType d100 = igh::d100<CalculationType>();
 143
 144                         int z;
 145                         if (lp_lat >=  d4044118) {          // 1|2
 146                           z = (lp_lon <= -d40 ? 1: 2);
 147                         }
 148                         else if (lp_lat >=  0) {            // 3|4
 149                           z = (lp_lon <= -d40 ? 3: 4);
 150                         }
 151                         else if (lp_lat >= -d4044118) {     // 5|6|7|8
 152                                if (lp_lon <= -d100) z =  5; // 5
 153                           else if (lp_lon <=  -d20) z =  6; // 6
 154                           else if (lp_lon <=   d80) z =  7; // 7
 155                           else z = 8;                       // 8
 156                         }
 157                         else {                              // 9|10|11|12
 158                                if (lp_lon <= -d100) z =  9; // 9
 159                           else if (lp_lon <=  -d20) z = 10; // 10
 160                           else if (lp_lon <=   d80) z = 11; // 11
 161                           else z = 12;                      // 12
 162                         }
 163
 164                         lp_lon -= this->m_proj_parm.pj[z-1]->params().lam0;
 165                         this->m_proj_parm.pj[z-1]->fwd(lp_lon, lp_lat, xy_x, xy_y);
 166                         xy_x += this->m_proj_parm.pj[z-1]->params().x0;
 167                         xy_y += this->m_proj_parm.pj[z-1]->params().y0;
 168                 }
 169
 170                 // INVERSE(s_inverse)  spheroid
 171                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 172                 inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
 173                 {
 174                     static const CalculationType d4044118 = igh::d4044118<CalculationType>();
 175                     static const CalculationType d10  =  igh::d10<CalculationType>();
 176                     static const CalculationType d20  =  igh::d20<CalculationType>();
 177                     static const CalculationType d40  =  igh::d40<CalculationType>();
 178                     static const CalculationType d50  =  igh::d50<CalculationType>();
 179                     static const CalculationType d60  =  igh::d60<CalculationType>();
 180                     static const CalculationType d80  =  igh::d80<CalculationType>();
 181                     static const CalculationType d90  =  igh::d90<CalculationType>();
 182                     static const CalculationType d100 = igh::d100<CalculationType>();
 183                     static const CalculationType d160 = igh::d160<CalculationType>();
 184                     static const CalculationType d180 = igh::d180<CalculationType>();
 185
 186                     static const CalculationType c2 = 2.0;
 187
 188                     const CalculationType y90 = this->m_proj_parm.dy0 + sqrt(c2); // lt=90 corresponds to y=y0+sqrt(2.0)
 189
 190                         int z = 0;
 191                         if (xy_y > y90+EPSLN || xy_y < -y90+EPSLN) // 0
 192                           z = 0;
 193                         else if (xy_y >=  d4044118)       // 1|2
 194                           z = (xy_x <= -d40? 1: 2);
 195                         else if (xy_y >=  0)              // 3|4
 196                           z = (xy_x <= -d40? 3: 4);
 197                         else if (xy_y >= -d4044118) {     // 5|6|7|8
 198                                if (xy_x <= -d100) z =  5; // 5
 199                           else if (xy_x <=  -d20) z =  6; // 6
 200                           else if (xy_x <=   d80) z =  7; // 7
 201                           else z = 8;                     // 8
 202                         }
 203                         else {                            // 9|10|11|12
 204                                if (xy_x <= -d100) z =  9; // 9
 205                           else if (xy_x <=  -d20) z = 10; // 10
 206                           else if (xy_x <=   d80) z = 11; // 11
 207                           else z = 12;                    // 12
 208                         }
 209
 210                         if (z)
 211                         {
 212                           int ok = 0;
 213
 214                           xy_x -= this->m_proj_parm.pj[z-1]->params().x0;
 215                           xy_y -= this->m_proj_parm.pj[z-1]->params().y0;
 216                           this->m_proj_parm.pj[z-1]->inv(xy_x, xy_y, lp_lon, lp_lat);
 217                           lp_lon += this->m_proj_parm.pj[z-1]->params().lam0;
 218
 219                           switch (z) {
 220                             case  1: ok = (lp_lon >= -d180-EPSLN && lp_lon <=  -d40+EPSLN) ||
 221                                          ((lp_lon >=  -d40-EPSLN && lp_lon <=  -d10+EPSLN) &&
 222                                           (lp_lat >=   d60-EPSLN && lp_lat <=   d90+EPSLN)); break;
 223                             case  2: ok = (lp_lon >=  -d40-EPSLN && lp_lon <=  d180+EPSLN) ||
 224                                          ((lp_lon >= -d180-EPSLN && lp_lon <= -d160+EPSLN) &&
 225                                           (lp_lat >=   d50-EPSLN && lp_lat <=   d90+EPSLN)) ||
 226                                          ((lp_lon >=  -d50-EPSLN && lp_lon <=  -d40+EPSLN) &&
 227                                           (lp_lat >=   d60-EPSLN && lp_lat <=   d90+EPSLN)); break;
 228                             case  3: ok = (lp_lon >= -d180-EPSLN && lp_lon <=  -d40+EPSLN); break;
 229                             case  4: ok = (lp_lon >=  -d40-EPSLN && lp_lon <=  d180+EPSLN); break;
 230                             case  5: ok = (lp_lon >= -d180-EPSLN && lp_lon <= -d100+EPSLN); break;
 231                             case  6: ok = (lp_lon >= -d100-EPSLN && lp_lon <=  -d20+EPSLN); break;
 232                             case  7: ok = (lp_lon >=  -d20-EPSLN && lp_lon <=   d80+EPSLN); break;
 233                             case  8: ok = (lp_lon >=   d80-EPSLN && lp_lon <=  d180+EPSLN); break;
 234                             case  9: ok = (lp_lon >= -d180-EPSLN && lp_lon <= -d100+EPSLN); break;
 235                             case 10: ok = (lp_lon >= -d100-EPSLN && lp_lon <=  -d20+EPSLN); break;
 236                             case 11: ok = (lp_lon >=  -d20-EPSLN && lp_lon <=   d80+EPSLN); break;
 237                             case 12: ok = (lp_lon >=   d80-EPSLN && lp_lon <=  d180+EPSLN); break;
 238                           }
 239
 240                           z = (!ok? 0: z); // projectable?
 241                         }
 242                      // if (!z) pj_errno = -15; // invalid x or y
 243                         if (!z) lp_lon = HUGE_VAL;
 244                         if (!z) lp_lat = HUGE_VAL;
 245                 }
 246
 247                 static inline std::string get_name()
 248                 {
 249                     return "igh_spheroid";
 250                 }
 251
 252             };
 253
 254             // Interrupted Goode Homolosine
 255             template <typename CalculationType, typename Parameters>
 256             inline void setup_igh(Parameters& par, par_igh<CalculationType, Parameters>& proj_parm)
 257             {
 258                 static const CalculationType d0   =  0;
 259                 static const CalculationType d4044118 = igh::d4044118<CalculationType>();
 260                 static const CalculationType d20  =  igh::d20<CalculationType>();
 261                 static const CalculationType d30  =  igh::d30<CalculationType>();
 262                 static const CalculationType d60  =  igh::d60<CalculationType>();
 263                 static const CalculationType d100 = igh::d100<CalculationType>();
 264                 static const CalculationType d140 = igh::d140<CalculationType>();
 265                 static const CalculationType d160 = igh::d160<CalculationType>();
 266
 267             /*
 268               Zones:
 269
 270                 -180            -40                       180
 271                   +--------------+-------------------------+    Zones 1,2,9,10,11 & 12:
 272                   |1             |2                        |      Mollweide projection
 273                   |              |                         |
 274                   +--------------+-------------------------+    Zones 3,4,5,6,7 & 8:
 275                   |3             |4                        |      Sinusoidal projection
 276                   |              |                         |
 277                 0 +-------+------+-+-----------+-----------+
 278                   |5      |6       |7          |8          |
 279                   |       |        |           |           |
 280                   +-------+--------+-----------+-----------+
 281                   |9      |10      |11         |12         |
 282                   |       |        |           |           |
 283                   +-------+--------+-----------+-----------+
 284                 -180    -100      -20         80          180
 285             */
 286
 287
 288                     CalculationType lp_lam = 0, lp_phi = d4044118;
 289                     CalculationType xy1_x, xy1_y;
 290                     CalculationType xy3_x, xy3_y;
 291
 292                     // sinusoidal zones
 293                     do_setup<sinu_entry>(3, par, proj_parm, -d100, d0, -d100);
 294                     do_setup<sinu_entry>(4, par, proj_parm,   d30, d0,   d30);
 295                     do_setup<sinu_entry>(5, par, proj_parm, -d160, d0, -d160);
 296                     do_setup<sinu_entry>(6, par, proj_parm,  -d60, d0,  -d60);
 297                     do_setup<sinu_entry>(7, par, proj_parm,   d20, d0,   d20);
 298                     do_setup<sinu_entry>(8, par, proj_parm,  d140, d0,  d140);
 299
 300                     // mollweide zones
 301                     do_setup<moll_entry>(1, 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, par, proj_parm,   d30,  proj_parm.dy0,   d30);
 313                     do_setup<moll_entry>( 9, par, proj_parm, -d160, -proj_parm.dy0, -d160);
 314                     do_setup<moll_entry>(10, par, proj_parm,  -d60, -proj_parm.dy0,  -d60);
 315                     do_setup<moll_entry>(11, par, proj_parm,   d20, -proj_parm.dy0,   d20);
 316                     do_setup<moll_entry>(12, par, proj_parm,  d140, -proj_parm.dy0,  d140);
 317
 318                     par.es = 0.;
 319             }
 320
 321     }} // namespace detail::igh
 322     #endif // doxygen
 323
 324     /*!
 325         \brief Interrupted Goode Homolosine projection
 326         \ingroup projections
 327         \tparam Geographic latlong point type
 328         \tparam Cartesian xy point type
 329         \tparam Parameters parameter type
 330         \par Projection characteristics
 331          - Pseudocylindrical
 332          - Spheroid
 333         \par Example
 334         \image html ex_igh.gif
 335     */
 336     template <typename CalculationType, typename Parameters>
 337     struct igh_spheroid : public detail::igh::base_igh_spheroid<CalculationType, Parameters>
 338     {
 339         inline igh_spheroid(const Parameters& par) : detail::igh::base_igh_spheroid<CalculationType, Parameters>(par)
 340         {
 341             detail::igh::setup_igh(this->m_par, this->m_proj_parm);
 342         }
 343     };
 344
 345     #ifndef DOXYGEN_NO_DETAIL
 346     namespace detail
 347     {
 348
 349         // Static projection
 350         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::igh, igh_spheroid, igh_spheroid)
 351
 352         // Factory entry(s)
 353         template <typename CalculationType, typename Parameters>
 354         class igh_entry : public detail::factory_entry<CalculationType, Parameters>
 355         {
 356             public :
 357                 virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const
 358                 {
 359                     return new base_v_fi<igh_spheroid<CalculationType, Parameters>, CalculationType, Parameters>(par);
 360                 }
 361         };
 362
 363         template <typename CalculationType, typename Parameters>
 364         inline void igh_init(detail::base_factory<CalculationType, Parameters>& factory)
 365         {
 366             factory.add_to_factory("igh", new igh_entry<CalculationType, Parameters>);
 367         }
 368
 369     } // namespace detail
 370     #endif // doxygen
 371
 372 } // namespace projections
 373
 374 }} // namespace boost::geometry
 375
 376 #endif // BOOST_GEOMETRY_PROJECTIONS_IGH_HPP
 377