boost/geometry/srs/projections/proj/gn_sinu.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_GN_SINU_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_GN_SINU_HPP
  42
  43 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
  44 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  45 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  46 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  47 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
  48 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  49 #include <boost/geometry/srs/projections/impl/projects.hpp>
  50
  51 #include <boost/geometry/util/math.hpp>
  52
  53 namespace boost { namespace geometry
  54 {
  55
  56 namespace projections
  57 {
  58     #ifndef DOXYGEN_NO_DETAIL
  59     namespace detail { namespace gn_sinu
  60     {
  61
  62             static const double epsilon10 = 1e-10;
  63             static const int max_iter = 8;
  64             static const double loop_tol = 1e-7;
  65
  66             template <typename T>
  67             struct par_gn_sinu_e
  68             {
  69                 detail::en<T> en;
  70             };
  71
  72             template <typename T>
  73             struct par_gn_sinu_s
  74             {
  75                 T m, n, C_x, C_y;
  76             };
  77
  78             /* Ellipsoidal Sinusoidal only */
  79
  80             template <typename T, typename Parameters>
  81             struct base_gn_sinu_ellipsoid
  82             {
  83                 par_gn_sinu_e<T> m_proj_parm;
  84
  85                 // FORWARD(e_forward)  ellipsoid
  86                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  87                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  88                 {
  89                     T s, c;
  90
  91                     xy_y = pj_mlfn(lp_lat, s = sin(lp_lat), c = cos(lp_lat), this->m_proj_parm.en);
  92                     xy_x = lp_lon * c / sqrt(1. - par.es * s * s);
  93                 }
  94
  95                 // INVERSE(e_inverse)  ellipsoid
  96                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  97                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
  98                 {
  99                     static const T half_pi = detail::half_pi<T>();
 100
 101                     T s;
 102
 103                     if ((s = fabs(lp_lat = pj_inv_mlfn(xy_y, par.es, this->m_proj_parm.en))) < half_pi) {
 104                         s = sin(lp_lat);
 105                         lp_lon = xy_x * sqrt(1. - par.es * s * s) / cos(lp_lat);
 106                     } else if ((s - epsilon10) < half_pi)
 107                         lp_lon = 0.;
 108                     else
 109                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 110                 }
 111                 /* General spherical sinusoidals */
 112
 113                 static inline std::string get_name()
 114                 {
 115                     return "gn_sinu_ellipsoid";
 116                 }
 117
 118             };
 119
 120             template <typename T, typename Parameters>
 121             struct base_gn_sinu_spheroid
 122             {
 123                 par_gn_sinu_s<T> m_proj_parm;
 124
 125                 // FORWARD(s_forward)  sphere
 126                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 127                 inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
 128                 {
 129                     if (this->m_proj_parm.m == 0.0)
 130                         lp_lat = this->m_proj_parm.n != 1. ? aasin(this->m_proj_parm.n * sin(lp_lat)): lp_lat;
 131                     else {
 132                         T k, V;
 133                         int i;
 134
 135                         k = this->m_proj_parm.n * sin(lp_lat);
 136                         for (i = max_iter; i ; --i) {
 137                             lp_lat -= V = (this->m_proj_parm.m * lp_lat + sin(lp_lat) - k) /
 138                                 (this->m_proj_parm.m + cos(lp_lat));
 139                             if (fabs(V) < loop_tol)
 140                                 break;
 141                         }
 142                         if (!i) {
 143                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 144                         }
 145                     }
 146                     xy_x = this->m_proj_parm.C_x * lp_lon * (this->m_proj_parm.m + cos(lp_lat));
 147                     xy_y = this->m_proj_parm.C_y * lp_lat;
 148                 }
 149
 150                 // INVERSE(s_inverse)  sphere
 151                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 152                 inline void inv(Parameters const& , T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 153                 {
 154                     xy_y /= this->m_proj_parm.C_y;
 155                     lp_lat = (this->m_proj_parm.m != 0.0) ? aasin((this->m_proj_parm.m * xy_y + sin(xy_y)) / this->m_proj_parm.n) :
 156                         ( this->m_proj_parm.n != 1. ? aasin(sin(xy_y) / this->m_proj_parm.n) : xy_y );
 157                     lp_lon = xy_x / (this->m_proj_parm.C_x * (this->m_proj_parm.m + cos(xy_y)));
 158                 }
 159
 160                 static inline std::string get_name()
 161                 {
 162                     return "gn_sinu_spheroid";
 163                 }
 164
 165             };
 166
 167             template <typename Parameters, typename T>
 168             inline void setup(Parameters& par, par_gn_sinu_s<T>& proj_parm)
 169             {
 170                 par.es = 0;
 171
 172                 proj_parm.C_x = (proj_parm.C_y = sqrt((proj_parm.m + 1.) / proj_parm.n))/(proj_parm.m + 1.);
 173             }
 174
 175
 176             // General Sinusoidal Series
 177             template <typename Params, typename Parameters, typename T>
 178             inline void setup_gn_sinu(Params const& params, Parameters& par, par_gn_sinu_s<T>& proj_parm)
 179             {
 180                 if (pj_param_f<srs::spar::n>(params, "n", srs::dpar::n, proj_parm.n)
 181                  && pj_param_f<srs::spar::m>(params, "m", srs::dpar::m, proj_parm.m)) {
 182                     if (proj_parm.n <= 0 || proj_parm.m < 0)
 183                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_m_or_n) );
 184                 } else
 185                     BOOST_THROW_EXCEPTION( projection_exception(error_invalid_m_or_n) );
 186
 187                 setup(par, proj_parm);
 188             }
 189
 190             // Sinusoidal (Sanson-Flamsteed)
 191             template <typename Parameters, typename T>
 192             inline void setup_sinu(Parameters const& par, par_gn_sinu_e<T>& proj_parm)
 193             {
 194                 proj_parm.en = pj_enfn<T>(par.es);
 195             }
 196
 197             // Sinusoidal (Sanson-Flamsteed)
 198             template <typename Parameters, typename T>
 199             inline void setup_sinu(Parameters& par, par_gn_sinu_s<T>& proj_parm)
 200             {
 201                 proj_parm.n = 1.;
 202                 proj_parm.m = 0.;
 203                 setup(par, proj_parm);
 204             }
 205
 206             // Eckert VI
 207             template <typename Parameters, typename T>
 208             inline void setup_eck6(Parameters& par, par_gn_sinu_s<T>& proj_parm)
 209             {
 210                 proj_parm.m = 1.;
 211                 proj_parm.n = 2.570796326794896619231321691;
 212                 setup(par, proj_parm);
 213             }
 214
 215             // McBryde-Thomas Flat-Polar Sinusoidal
 216             template <typename Parameters, typename T>
 217             inline void setup_mbtfps(Parameters& par, par_gn_sinu_s<T>& proj_parm)
 218             {
 219                 proj_parm.m = 0.5;
 220                 proj_parm.n = 1.785398163397448309615660845;
 221                 setup(par, proj_parm);
 222             }
 223
 224     }} // namespace detail::gn_sinu
 225     #endif // doxygen
 226
 227     /*!
 228         \brief General Sinusoidal Series projection
 229         \ingroup projections
 230         \tparam Geographic latlong point type
 231         \tparam Cartesian xy point type
 232         \tparam Parameters parameter type
 233         \par Projection characteristics
 234          - Pseudocylindrical
 235          - Spheroid
 236         \par Projection parameters
 237          - m (real)
 238          - n (real)
 239         \par Example
 240         \image html ex_gn_sinu.gif
 241     */
 242     template <typename T, typename Parameters>
 243     struct gn_sinu_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
 244     {
 245         template <typename Params>
 246         inline gn_sinu_spheroid(Params const& params, Parameters & par)
 247         {
 248             detail::gn_sinu::setup_gn_sinu(params, par, this->m_proj_parm);
 249         }
 250     };
 251
 252     /*!
 253         \brief Sinusoidal (Sanson-Flamsteed) projection
 254         \ingroup projections
 255         \tparam Geographic latlong point type
 256         \tparam Cartesian xy point type
 257         \tparam Parameters parameter type
 258         \par Projection characteristics
 259          - Pseudocylindrical
 260          - Spheroid
 261          - Ellipsoid
 262         \par Example
 263         \image html ex_sinu.gif
 264     */
 265     template <typename T, typename Parameters>
 266     struct sinu_ellipsoid : public detail::gn_sinu::base_gn_sinu_ellipsoid<T, Parameters>
 267     {
 268         template <typename Params>
 269         inline sinu_ellipsoid(Params const& , Parameters & par)
 270         {
 271             detail::gn_sinu::setup_sinu(par, this->m_proj_parm);
 272         }
 273     };
 274
 275     /*!
 276         \brief Sinusoidal (Sanson-Flamsteed) projection
 277         \ingroup projections
 278         \tparam Geographic latlong point type
 279         \tparam Cartesian xy point type
 280         \tparam Parameters parameter type
 281         \par Projection characteristics
 282          - Pseudocylindrical
 283          - Spheroid
 284          - Ellipsoid
 285         \par Example
 286         \image html ex_sinu.gif
 287     */
 288     template <typename T, typename Parameters>
 289     struct sinu_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
 290     {
 291         template <typename Params>
 292         inline sinu_spheroid(Params const& , Parameters & par)
 293         {
 294             detail::gn_sinu::setup_sinu(par, this->m_proj_parm);
 295         }
 296     };
 297
 298     /*!
 299         \brief Eckert VI projection
 300         \ingroup projections
 301         \tparam Geographic latlong point type
 302         \tparam Cartesian xy point type
 303         \tparam Parameters parameter type
 304         \par Projection characteristics
 305          - Pseudocylindrical
 306          - Spheroid
 307         \par Example
 308         \image html ex_eck6.gif
 309     */
 310     template <typename T, typename Parameters>
 311     struct eck6_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
 312     {
 313         template <typename Params>
 314         inline eck6_spheroid(Params const& , Parameters & par)
 315         {
 316             detail::gn_sinu::setup_eck6(par, this->m_proj_parm);
 317         }
 318     };
 319
 320     /*!
 321         \brief McBryde-Thomas Flat-Polar Sinusoidal projection
 322         \ingroup projections
 323         \tparam Geographic latlong point type
 324         \tparam Cartesian xy point type
 325         \tparam Parameters parameter type
 326         \par Projection characteristics
 327          - Pseudocylindrical
 328          - Spheroid
 329         \par Example
 330         \image html ex_mbtfps.gif
 331     */
 332     template <typename T, typename Parameters>
 333     struct mbtfps_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
 334     {
 335         template <typename Params>
 336         inline mbtfps_spheroid(Params const& , Parameters & par)
 337         {
 338             detail::gn_sinu::setup_mbtfps(par, this->m_proj_parm);
 339         }
 340     };
 341
 342     #ifndef DOXYGEN_NO_DETAIL
 343     namespace detail
 344     {
 345
 346         // Static projection
 347         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_gn_sinu, gn_sinu_spheroid)
 348         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_sinu, sinu_spheroid, sinu_ellipsoid)
 349         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_eck6, eck6_spheroid)
 350         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_mbtfps, mbtfps_spheroid)
 351
 352         // Factory entry(s)
 353         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(gn_sinu_entry, gn_sinu_spheroid)
 354         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(sinu_entry, sinu_spheroid, sinu_ellipsoid)
 355         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(eck6_entry, eck6_spheroid)
 356         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(mbtfps_entry, mbtfps_spheroid)
 357
 358         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(gn_sinu_init)
 359         {
 360             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(gn_sinu, gn_sinu_entry);
 361             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(sinu, sinu_entry);
 362             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(eck6, eck6_entry);
 363             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(mbtfps, mbtfps_entry);
 364         }
 365
 366     } // namespace detail
 367     #endif // doxygen
 368
 369 } // namespace projections
 370
 371 }} // namespace boost::geometry
 372
 373 #endif // BOOST_GEOMETRY_PROJECTIONS_GN_SINU_HPP
 374