projections/proj/igh.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0005
0006 // This file was modified by Oracle on 2017-2021.
0007 // Modifications copyright (c) 2017-2021, Oracle and/or its affiliates.
0008 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0009
0010 // Use, modification and distribution is subject to the Boost Software License,
0011 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0012 // http://www.boost.org/LICENSE_1_0.txt)
0013
0014 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0015 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0016 // PROJ4 is maintained by Frank Warmerdam
0017 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0018
0019 // Last updated version of proj: 5.0.0
0020
0021 // Original copyright notice:
0022
0023 // Permission is hereby granted, free of charge, to any person obtaining a
0024 // copy of this software and associated documentation files (the "Software"),
0025 // to deal in the Software without restriction, including without limitation
0026 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0027 // and/or sell copies of the Software, and to permit persons to whom the
0028 // Software is furnished to do so, subject to the following conditions:
0029
0030 // The above copyright notice and this permission notice shall be included
0031 // in all copies or substantial portions of the Software.
0032
0033 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0034 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0035 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0036 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0037 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0038 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0039 // DEALINGS IN THE SOFTWARE.
0040
0041 #ifndef BOOST_GEOMETRY_PROJECTIONS_IGH_HPP
0042 #define BOOST_GEOMETRY_PROJECTIONS_IGH_HPP
0043
0044 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0045 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0046 #include <boost/geometry/srs/projections/impl/projects.hpp>
0047 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0048 #include <boost/geometry/srs/projections/proj/gn_sinu.hpp>
0049 #include <boost/geometry/srs/projections/proj/moll.hpp>
0050
0051 #include <boost/geometry/util/math.hpp>
0052
0053 namespace boost { namespace geometry
0054 {
0055
0056 namespace projections
0057 {
0058     #ifndef DOXYGEN_NO_DETAIL
0059     namespace detail { namespace igh
0060     {
0061             template <typename T>
0062             struct par_igh_zone
0063             {
0064                 T x0;
0065                 T y0;
0066                 T lam0;
0067             };
0068
0069             // NOTE: x0, y0, lam0 are not used in moll nor sinu projections
0070             // so it is a waste of memory to keep 12 copies of projections
0071             // with parameters as in the original Proj4.
0072
0073             // TODO: It would be possible to further decrease the size of par_igh
0074             // because spherical sinu and moll has constant parameters.
0075             // TODO: Furthermore there is no need to store par_igh_zone parameters
0076             // since they are constant for zones. In both fwd() and inv() there are
0077             // parts of code dependent on specific zones (if statements) anyway
0078             // so these parameters could be hardcoded there instead of stored.
0079
0080             template <typename T, typename Parameters>
0081             struct par_igh
0082             {
0083                 moll_spheroid<T, Parameters> moll;
0084                 sinu_spheroid<T, Parameters> sinu;
0085                 par_igh_zone<T> zones[12];
0086                 T dy0;
0087
0088                 // NOTE: The constructors of moll and sinu projections sets
0089                 // par.es = 0
0090
0091                 template <typename Params>
0092                 inline par_igh(Params const& params, Parameters & par)
0093                     : moll(params, par)
0094                     , sinu(params, par)
0095                 {}
0096
0097                 inline void fwd(int zone, Parameters const& par, T const& lp_lon, T const& lp_lat, T & xy_x, T & xy_y) const
0098                 {
0099                     if (zone <= 2 || zone >= 9) // 1, 2, 9, 10, 11, 12
0100                         moll.fwd(par, lp_lon, lp_lat, xy_x, xy_y);
0101                     else // 3, 4, 5, 6, 7, 8
0102                         sinu.fwd(par, lp_lon, lp_lat, xy_x, xy_y);
0103                 }
0104
0105                 inline void inv(int zone, Parameters const& par, T const& xy_x, T const& xy_y, T & lp_lon, T & lp_lat) const
0106                 {
0107                     if (zone <= 2 || zone >= 9) // 1, 2, 9, 10, 11, 12
0108                         moll.inv(par, xy_x, xy_y, lp_lon, lp_lat);
0109                     else // 3, 4, 5, 6, 7, 8
0110                         sinu.inv(par, xy_x, xy_y, lp_lon, lp_lat);
0111                 }
0112
0113                 inline void set_zone(int zone, T const& x_0, T const& y_0, T const& lon_0)
0114                 {
0115                     zones[zone - 1].x0 = x_0;
0116                     zones[zone - 1].y0 = y_0;
0117                     zones[zone - 1].lam0 = lon_0;
0118                 }
0119
0120                 inline par_igh_zone<T> const& get_zone(int zone) const
0121                 {
0122                     return zones[zone - 1];
0123                 }
0124             };
0125
0126             /* 40d 44' 11.8" [degrees] */
0127             template <typename T>
0128             inline T d4044118() { return (T(40) + T(44)/T(60.) + T(11.8)/T(3600.)) * geometry::math::d2r<T>(); }
0129
0130             template <typename T>
0131             inline T d10() { return T(10) * geometry::math::d2r<T>(); }
0132             template <typename T>
0133             inline T d20() { return T(20) * geometry::math::d2r<T>(); }
0134             template <typename T>
0135             inline T d30() { return T(30) * geometry::math::d2r<T>(); }
0136             template <typename T>
0137             inline T d40() { return T(40) * geometry::math::d2r<T>(); }
0138             template <typename T>
0139             inline T d50() { return T(50) * geometry::math::d2r<T>(); }
0140             template <typename T>
0141             inline T d60() { return T(60) * geometry::math::d2r<T>(); }
0142             template <typename T>
0143             inline T d80() { return T(80) * geometry::math::d2r<T>(); }
0144             template <typename T>
0145             inline T d90() { return T(90) * geometry::math::d2r<T>(); }
0146             template <typename T>
0147             inline T d100() { return T(100) * geometry::math::d2r<T>(); }
0148             template <typename T>
0149             inline T d140() { return T(140) * geometry::math::d2r<T>(); }
0150             template <typename T>
0151             inline T d160() { return T(160) * geometry::math::d2r<T>(); }
0152             template <typename T>
0153             inline T d180() { return T(180) * geometry::math::d2r<T>(); }
0154
0155             static const double epsilon = 1.e-10; // allow a little 'slack' on zone edge positions
0156
0157             template <typename T, typename Parameters>
0158             struct base_igh_spheroid
0159             {
0160                 par_igh<T, Parameters> m_proj_parm;
0161
0162                 template <typename Params>
0163                 inline base_igh_spheroid(Params const& params, Parameters & par)
0164                     : m_proj_parm(params, par)
0165                 {}
0166
0167                 // FORWARD(s_forward)  spheroid
0168                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0169                 inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0170                 {
0171                     static const T d4044118 = igh::d4044118<T>();
0172                     static const T d20  =  igh::d20<T>();
0173                     static const T d40  =  igh::d40<T>();
0174                     static const T d80  =  igh::d80<T>();
0175                     static const T d100 = igh::d100<T>();
0176
0177                         int z;
0178                         if (lp_lat >=  d4044118) {          // 1|2
0179                           z = (lp_lon <= -d40 ? 1: 2);
0180                         }
0181                         else if (lp_lat >=  0) {            // 3|4
0182                           z = (lp_lon <= -d40 ? 3: 4);
0183                         }
0184                         else if (lp_lat >= -d4044118) {     // 5|6|7|8
0185                                if (lp_lon <= -d100) z =  5; // 5
0186                           else if (lp_lon <=  -d20) z =  6; // 6
0187                           else if (lp_lon <=   d80) z =  7; // 7
0188                           else z = 8;                       // 8
0189                         }
0190                         else {                              // 9|10|11|12
0191                                if (lp_lon <= -d100) z =  9; // 9
0192                           else if (lp_lon <=  -d20) z = 10; // 10
0193                           else if (lp_lon <=   d80) z = 11; // 11
0194                           else z = 12;                      // 12
0195                         }
0196
0197                         lp_lon -= this->m_proj_parm.get_zone(z).lam0;
0198                         this->m_proj_parm.fwd(z, par, lp_lon, lp_lat, xy_x, xy_y);
0199                         xy_x += this->m_proj_parm.get_zone(z).x0;
0200                         xy_y += this->m_proj_parm.get_zone(z).y0;
0201                 }
0202
0203                 // INVERSE(s_inverse)  spheroid
0204                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0205                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0206                 {
0207                     static const T d4044118 = igh::d4044118<T>();
0208                     static const T d10  =  igh::d10<T>();
0209                     static const T d20  =  igh::d20<T>();
0210                     static const T d40  =  igh::d40<T>();
0211                     static const T d50  =  igh::d50<T>();
0212                     static const T d60  =  igh::d60<T>();
0213                     static const T d80  =  igh::d80<T>();
0214                     static const T d90  =  igh::d90<T>();
0215                     static const T d100 = igh::d100<T>();
0216                     static const T d160 = igh::d160<T>();
0217                     static const T d180 = igh::d180<T>();
0218
0219                     static const T c2 = 2.0;
0220
0221                     const T y90 = this->m_proj_parm.dy0 + sqrt(c2); // lt=90 corresponds to y=y0+sqrt(2.0)
0222
0223                         int z = 0;
0224                         if (xy_y > y90+epsilon || xy_y < -y90+epsilon) // 0
0225                           z = 0;
0226                         else if (xy_y >=  d4044118)       // 1|2
0227                           z = (xy_x <= -d40? 1: 2);
0228                         else if (xy_y >=  0)              // 3|4
0229                           z = (xy_x <= -d40? 3: 4);
0230                         else if (xy_y >= -d4044118) {     // 5|6|7|8
0231                                if (xy_x <= -d100) z =  5; // 5
0232                           else if (xy_x <=  -d20) z =  6; // 6
0233                           else if (xy_x <=   d80) z =  7; // 7
0234                           else z = 8;                     // 8
0235                         }
0236                         else {                            // 9|10|11|12
0237                                if (xy_x <= -d100) z =  9; // 9
0238                           else if (xy_x <=  -d20) z = 10; // 10
0239                           else if (xy_x <=   d80) z = 11; // 11
0240                           else z = 12;                    // 12
0241                         }
0242
0243                         if (z)
0244                         {
0245                           int ok = 0;
0246
0247                           xy_x -= this->m_proj_parm.get_zone(z).x0;
0248                           xy_y -= this->m_proj_parm.get_zone(z).y0;
0249                           this->m_proj_parm.inv(z, par, xy_x, xy_y, lp_lon, lp_lat);
0250                           lp_lon += this->m_proj_parm.get_zone(z).lam0;
0251
0252                           switch (z) {
0253                             case  1: ok = (lp_lon >= -d180-epsilon && lp_lon <=  -d40+epsilon) ||
0254                                          ((lp_lon >=  -d40-epsilon && lp_lon <=  -d10+epsilon) &&
0255                                           (lp_lat >=   d60-epsilon && lp_lat <=   d90+epsilon)); break;
0256                             case  2: ok = (lp_lon >=  -d40-epsilon && lp_lon <=  d180+epsilon) ||
0257                                          ((lp_lon >= -d180-epsilon && lp_lon <= -d160+epsilon) &&
0258                                           (lp_lat >=   d50-epsilon && lp_lat <=   d90+epsilon)) ||
0259                                          ((lp_lon >=  -d50-epsilon && lp_lon <=  -d40+epsilon) &&
0260                                           (lp_lat >=   d60-epsilon && lp_lat <=   d90+epsilon)); break;
0261                             case  3: ok = (lp_lon >= -d180-epsilon && lp_lon <=  -d40+epsilon); break;
0262                             case  4: ok = (lp_lon >=  -d40-epsilon && lp_lon <=  d180+epsilon); break;
0263                             case  5: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d100+epsilon); break;
0264                             case  6: ok = (lp_lon >= -d100-epsilon && lp_lon <=  -d20+epsilon); break;
0265                             case  7: ok = (lp_lon >=  -d20-epsilon && lp_lon <=   d80+epsilon); break;
0266                             case  8: ok = (lp_lon >=   d80-epsilon && lp_lon <=  d180+epsilon); break;
0267                             case  9: ok = (lp_lon >= -d180-epsilon && lp_lon <= -d100+epsilon); break;
0268                             case 10: ok = (lp_lon >= -d100-epsilon && lp_lon <=  -d20+epsilon); break;
0269                             case 11: ok = (lp_lon >=  -d20-epsilon && lp_lon <=   d80+epsilon); break;
0270                             case 12: ok = (lp_lon >=   d80-epsilon && lp_lon <=  d180+epsilon); break;
0271                           }
0272
0273                           z = (!ok? 0: z); // projectable?
0274                         }
0275                      // if (!z) pj_errno = -15; // invalid x or y
0276                         if (!z) lp_lon = HUGE_VAL;
0277                         if (!z) lp_lat = HUGE_VAL;
0278                 }
0279
0280                 static inline std::string get_name()
0281                 {
0282                     return "igh_spheroid";
0283                 }
0284
0285             };
0286
0287             // Interrupted Goode Homolosine
0288             template <typename Params, typename Parameters, typename T>
0289             inline void setup_igh(Params const& , Parameters& par, par_igh<T, Parameters>& proj_parm)
0290             {
0291                 static const T d0   =  0;
0292                 static const T d4044118 = igh::d4044118<T>();
0293                 static const T d20  =  igh::d20<T>();
0294                 static const T d30  =  igh::d30<T>();
0295                 static const T d60  =  igh::d60<T>();
0296                 static const T d100 = igh::d100<T>();
0297                 static const T d140 = igh::d140<T>();
0298                 static const T d160 = igh::d160<T>();
0299
0300             /*
0301               Zones:
0302
0303                 -180            -40                       180
0304                   +--------------+-------------------------+    Zones 1,2,9,10,11 & 12:
0305                   |1             |2                        |      Mollweide projection
0306                   |              |                         |
0307                   +--------------+-------------------------+    Zones 3,4,5,6,7 & 8:
0308                   |3             |4                        |      Sinusoidal projection
0309                   |              |                         |
0310                 0 +-------+------+-+-----------+-----------+
0311                   |5      |6       |7          |8          |
0312                   |       |        |           |           |
0313                   +-------+--------+-----------+-----------+
0314                   |9      |10      |11         |12         |
0315                   |       |        |           |           |
0316                   +-------+--------+-----------+-----------+
0317                 -180    -100      -20         80          180
0318             */
0319
0320                     T lp_lam = 0, lp_phi = d4044118;
0321                     T xy1_x, xy1_y;
0322                     T xy3_x, xy3_y;
0323
0324                     // sinusoidal zones
0325                     proj_parm.set_zone(3, -d100, d0, -d100);
0326                     proj_parm.set_zone(4,   d30, d0,   d30);
0327                     proj_parm.set_zone(5, -d160, d0, -d160);
0328                     proj_parm.set_zone(6,  -d60, d0,  -d60);
0329                     proj_parm.set_zone(7,   d20, d0,   d20);
0330                     proj_parm.set_zone(8,  d140, d0,  d140);
0331
0332                     // mollweide zones
0333                     proj_parm.set_zone(1, -d100, d0, -d100);
0334
0335                     // NOTE: x0, y0, lam0 are not used in moll nor sinu fwd
0336                     // so the order of initialization doesn't matter that much.
0337                     // But keep the original one from Proj4.
0338
0339                     // y0 ?
0340                     proj_parm.fwd(1, par, lp_lam, lp_phi, xy1_x, xy1_y); // zone 1
0341                     proj_parm.fwd(3, par, lp_lam, lp_phi, xy3_x, xy3_y); // zone 3
0342                     // y0 + xy1_y = xy3_y for lt = 40d44'11.8"
0343                     proj_parm.dy0 = xy3_y - xy1_y;
0344
0345                     proj_parm.zones[0].y0 = proj_parm.dy0; // zone 1
0346
0347                     // mollweide zones (cont'd)
0348                     proj_parm.set_zone(2,   d30,  proj_parm.dy0,   d30);
0349                     proj_parm.set_zone(9, -d160, -proj_parm.dy0, -d160);
0350                     proj_parm.set_zone(10, -d60, -proj_parm.dy0,  -d60);
0351                     proj_parm.set_zone(11,  d20, -proj_parm.dy0,   d20);
0352                     proj_parm.set_zone(12, d140, -proj_parm.dy0,  d140);
0353
0354                     // NOTE: Already done before in sinu and moll constructor
0355                     //par.es = 0.;
0356             }
0357
0358     }} // namespace detail::igh
0359     #endif // doxygen
0360
0361     /*!
0362         \brief Interrupted Goode Homolosine projection
0363         \ingroup projections
0364         \tparam Geographic latlong point type
0365         \tparam Cartesian xy point type
0366         \tparam Parameters parameter type
0367         \par Projection characteristics
0368          - Pseudocylindrical
0369          - Spheroid
0370         \par Example
0371         \image html ex_igh.gif
0372     */
0373     template <typename T, typename Parameters>
0374     struct igh_spheroid : public detail::igh::base_igh_spheroid<T, Parameters>
0375     {
0376         template <typename Params>
0377         inline igh_spheroid(Params const& params, Parameters & par)
0378             : detail::igh::base_igh_spheroid<T, Parameters>(params, par)
0379         {
0380             detail::igh::setup_igh(params, par, this->m_proj_parm);
0381         }
0382     };
0383
0384     #ifndef DOXYGEN_NO_DETAIL
0385     namespace detail
0386     {
0387
0388         // Static projection
0389         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_igh, igh_spheroid)
0390
0391         // Factory entry(s)
0392         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(igh_entry, igh_spheroid)
0393
0394         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(igh_init)
0395         {
0396             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(igh, igh_entry)
0397         }
0398
0399     } // namespace detail
0400     #endif // doxygen
0401
0402 } // namespace projections
0403
0404 }} // namespace boost::geometry
0405
0406 #endif // BOOST_GEOMETRY_PROJECTIONS_IGH_HPP
0407