boost/geometry/srs/projections/proj/sconics.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_SCONICS_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_SCONICS_HPP
  42
  43
  44 #include <boost/geometry/util/math.hpp>
  45 #include <boost/math/special_functions/hypot.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/factory_entry.hpp>
  50 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  51 #include <boost/geometry/srs/projections/impl/projects.hpp>
  52
  53 namespace boost { namespace geometry
  54 {
  55
  56 namespace projections
  57 {
  58     #ifndef DOXYGEN_NO_DETAIL
  59     namespace detail { namespace sconics
  60     {
  61
  62             enum proj_type {
  63                 proj_euler  = 0,
  64                 proj_murd1  = 1,
  65                 proj_murd2  = 2,
  66                 proj_murd3  = 3,
  67                 proj_pconic = 4,
  68                 proj_tissot = 5,
  69                 proj_vitk1  = 6
  70             };
  71             static const double epsilon10 = 1.e-10;
  72             static const double epsilon = 1e-10;
  73
  74             template <typename T>
  75             struct par_sconics
  76             {
  77                 T   n;
  78                 T   rho_c;
  79                 T   rho_0;
  80                 T   sig;
  81                 T   c1, c2;
  82                 proj_type type;
  83             };
  84
  85             /* get common factors for simple conics */
  86             template <typename Params, typename T>
  87             inline int phi12(Params const& params, par_sconics<T>& proj_parm, T *del)
  88             {
  89                 T p1, p2;
  90                 int err = 0;
  91
  92                 if (!pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, p1) ||
  93                     !pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, p2)) {
  94                     err = -41;
  95                 } else {
  96                     //p1 = pj_get_param_r(par.params, "lat_1"); // set above
  97                     //p2 = pj_get_param_r(par.params, "lat_2"); // set above
  98                     *del = 0.5 * (p2 - p1);
  99                     proj_parm.sig = 0.5 * (p2 + p1);
 100                     err = (fabs(*del) < epsilon || fabs(proj_parm.sig) < epsilon) ? -42 : 0;
 101                 }
 102                 return err;
 103             }
 104
 105             template <typename T, typename Parameters>
 106             struct base_sconics_spheroid
 107             {
 108                 par_sconics<T> m_proj_parm;
 109
 110                 // FORWARD(s_forward)  spheroid
 111                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 112                 inline void fwd(Parameters const& , T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 113                 {
 114                     T rho;
 115
 116                     switch (this->m_proj_parm.type) {
 117                     case proj_murd2:
 118                         rho = this->m_proj_parm.rho_c + tan(this->m_proj_parm.sig - lp_lat);
 119                         break;
 120                     case proj_pconic:
 121                         rho = this->m_proj_parm.c2 * (this->m_proj_parm.c1 - tan(lp_lat - this->m_proj_parm.sig));
 122                         break;
 123                     default:
 124                         rho = this->m_proj_parm.rho_c - lp_lat;
 125                         break;
 126                     }
 127                     xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
 128                     xy_y = this->m_proj_parm.rho_0 - rho * cos(lp_lon);
 129                 }
 130
 131                 // INVERSE(s_inverse)  ellipsoid & spheroid
 132                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 133                 inline void inv(Parameters const& , T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 134                 {
 135                     T rho;
 136
 137                     rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho_0 - xy_y);
 138                     if (this->m_proj_parm.n < 0.) {
 139                         rho = - rho;
 140                         xy_x = - xy_x;
 141                         xy_y = - xy_y;
 142                     }
 143
 144                     lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
 145
 146                     switch (this->m_proj_parm.type) {
 147                     case proj_pconic:
 148                         lp_lat = atan(this->m_proj_parm.c1 - rho / this->m_proj_parm.c2) + this->m_proj_parm.sig;
 149                         break;
 150                     case proj_murd2:
 151                         lp_lat = this->m_proj_parm.sig - atan(rho - this->m_proj_parm.rho_c);
 152                         break;
 153                     default:
 154                         lp_lat = this->m_proj_parm.rho_c - rho;
 155                     }
 156                 }
 157
 158                 static inline std::string get_name()
 159                 {
 160                     return "sconics_spheroid";
 161                 }
 162
 163             };
 164
 165             template <typename Params, typename Parameters, typename T>
 166             inline void setup(Params const& params, Parameters& par, par_sconics<T>& proj_parm, proj_type type)
 167             {
 168                 static const T half_pi = detail::half_pi<T>();
 169
 170                 T del, cs;
 171                 int err;
 172
 173                 proj_parm.type = type;
 174
 175                 err = phi12(params, proj_parm, &del);
 176                 if(err)
 177                     BOOST_THROW_EXCEPTION( projection_exception(err) );
 178
 179                 switch (proj_parm.type) {
 180                 case proj_tissot:
 181                     proj_parm.n = sin(proj_parm.sig);
 182                     cs = cos(del);
 183                     proj_parm.rho_c = proj_parm.n / cs + cs / proj_parm.n;
 184                     proj_parm.rho_0 = sqrt((proj_parm.rho_c - 2 * sin(par.phi0))/proj_parm.n);
 185                     break;
 186                 case proj_murd1:
 187                     proj_parm.rho_c = sin(del)/(del * tan(proj_parm.sig)) + proj_parm.sig;
 188                     proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
 189                     proj_parm.n = sin(proj_parm.sig);
 190                     break;
 191                 case proj_murd2:
 192                     proj_parm.rho_c = (cs = sqrt(cos(del))) / tan(proj_parm.sig);
 193                     proj_parm.rho_0 = proj_parm.rho_c + tan(proj_parm.sig - par.phi0);
 194                     proj_parm.n = sin(proj_parm.sig) * cs;
 195                     break;
 196                 case proj_murd3:
 197                     proj_parm.rho_c = del / (tan(proj_parm.sig) * tan(del)) + proj_parm.sig;
 198                     proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
 199                     proj_parm.n = sin(proj_parm.sig) * sin(del) * tan(del) / (del * del);
 200                     break;
 201                 case proj_euler:
 202                     proj_parm.n = sin(proj_parm.sig) * sin(del) / del;
 203                     del *= 0.5;
 204                     proj_parm.rho_c = del / (tan(del) * tan(proj_parm.sig)) + proj_parm.sig;
 205                     proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
 206                     break;
 207                 case proj_pconic:
 208                     proj_parm.n = sin(proj_parm.sig);
 209                     proj_parm.c2 = cos(del);
 210                     proj_parm.c1 = 1./tan(proj_parm.sig);
 211                     if (fabs(del = par.phi0 - proj_parm.sig) - epsilon10 >= half_pi)
 212                         BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_half_pi_from_mean) );
 213                     proj_parm.rho_0 = proj_parm.c2 * (proj_parm.c1 - tan(del));
 214                     break;
 215                 case proj_vitk1:
 216                     proj_parm.n = (cs = tan(del)) * sin(proj_parm.sig) / del;
 217                     proj_parm.rho_c = del / (cs * tan(proj_parm.sig)) + proj_parm.sig;
 218                     proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
 219                     break;
 220                 }
 221
 222                 par.es = 0;
 223             }
 224
 225
 226             // Euler
 227             template <typename Params, typename Parameters, typename T>
 228             inline void setup_euler(Params const& params, Parameters& par, par_sconics<T>& proj_parm)
 229             {
 230                 setup(params, par, proj_parm, proj_euler);
 231             }
 232
 233             // Tissot
 234             template <typename Params, typename Parameters, typename T>
 235             inline void setup_tissot(Params const& params, Parameters& par, par_sconics<T>& proj_parm)
 236             {
 237                 setup(params, par, proj_parm, proj_tissot);
 238             }
 239
 240             // Murdoch I
 241             template <typename Params, typename Parameters, typename T>
 242             inline void setup_murd1(Params const& params, Parameters& par, par_sconics<T>& proj_parm)
 243             {
 244                 setup(params, par, proj_parm, proj_murd1);
 245             }
 246
 247             // Murdoch II
 248             template <typename Params, typename Parameters, typename T>
 249             inline void setup_murd2(Params const& params, Parameters& par, par_sconics<T>& proj_parm)
 250             {
 251                 setup(params, par, proj_parm, proj_murd2);
 252             }
 253
 254             // Murdoch III
 255             template <typename Params, typename Parameters, typename T>
 256             inline void setup_murd3(Params const& params, Parameters& par, par_sconics<T>& proj_parm)
 257             {
 258                 setup(params, par, proj_parm, proj_murd3);
 259             }
 260
 261             // Perspective Conic
 262             template <typename Params, typename Parameters, typename T>
 263             inline void setup_pconic(Params const& params, Parameters& par, par_sconics<T>& proj_parm)
 264             {
 265                 setup(params, par, proj_parm, proj_pconic);
 266             }
 267
 268             // Vitkovsky I
 269             template <typename Params, typename Parameters, typename T>
 270             inline void setup_vitk1(Params const& params, Parameters& par, par_sconics<T>& proj_parm)
 271             {
 272                 setup(params, par, proj_parm, proj_vitk1);
 273             }
 274
 275     }} // namespace detail::sconics
 276     #endif // doxygen
 277
 278     /*!
 279         \brief Tissot projection
 280         \ingroup projections
 281         \tparam Geographic latlong point type
 282         \tparam Cartesian xy point type
 283         \tparam Parameters parameter type
 284         \par Projection characteristics
 285          - Conic
 286          - Spheroid
 287         \par Projection parameters
 288          - lat_1: Latitude of first standard parallel
 289          - lat_2: Latitude of second standard parallel
 290         \par Example
 291         \image html ex_tissot.gif
 292     */
 293     template <typename T, typename Parameters>
 294     struct tissot_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
 295     {
 296         template <typename Params>
 297         inline tissot_spheroid(Params const& params, Parameters& par)
 298         {
 299             detail::sconics::setup_tissot(params, par, this->m_proj_parm);
 300         }
 301     };
 302
 303     /*!
 304         \brief Murdoch I projection
 305         \ingroup projections
 306         \tparam Geographic latlong point type
 307         \tparam Cartesian xy point type
 308         \tparam Parameters parameter type
 309         \par Projection characteristics
 310          - Conic
 311          - Spheroid
 312         \par Projection parameters
 313          - lat_1: Latitude of first standard parallel
 314          - lat_2: Latitude of second standard parallel
 315         \par Example
 316         \image html ex_murd1.gif
 317     */
 318     template <typename T, typename Parameters>
 319     struct murd1_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
 320     {
 321         template <typename Params>
 322         inline murd1_spheroid(Params const& params, Parameters& par)
 323         {
 324             detail::sconics::setup_murd1(params, par, this->m_proj_parm);
 325         }
 326     };
 327
 328     /*!
 329         \brief Murdoch II projection
 330         \ingroup projections
 331         \tparam Geographic latlong point type
 332         \tparam Cartesian xy point type
 333         \tparam Parameters parameter type
 334         \par Projection characteristics
 335          - Conic
 336          - Spheroid
 337         \par Projection parameters
 338          - lat_1: Latitude of first standard parallel
 339          - lat_2: Latitude of second standard parallel
 340         \par Example
 341         \image html ex_murd2.gif
 342     */
 343     template <typename T, typename Parameters>
 344     struct murd2_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
 345     {
 346         template <typename Params>
 347         inline murd2_spheroid(Params const& params, Parameters& par)
 348         {
 349             detail::sconics::setup_murd2(params, par, this->m_proj_parm);
 350         }
 351     };
 352
 353     /*!
 354         \brief Murdoch III projection
 355         \ingroup projections
 356         \tparam Geographic latlong point type
 357         \tparam Cartesian xy point type
 358         \tparam Parameters parameter type
 359         \par Projection characteristics
 360          - Conic
 361          - Spheroid
 362         \par Projection parameters
 363          - lat_1: Latitude of first standard parallel
 364          - lat_2: Latitude of second standard parallel
 365         \par Example
 366         \image html ex_murd3.gif
 367     */
 368     template <typename T, typename Parameters>
 369     struct murd3_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
 370     {
 371         template <typename Params>
 372         inline murd3_spheroid(Params const& params, Parameters& par)
 373         {
 374             detail::sconics::setup_murd3(params, par, this->m_proj_parm);
 375         }
 376     };
 377
 378     /*!
 379         \brief Euler projection
 380         \ingroup projections
 381         \tparam Geographic latlong point type
 382         \tparam Cartesian xy point type
 383         \tparam Parameters parameter type
 384         \par Projection characteristics
 385          - Conic
 386          - Spheroid
 387         \par Projection parameters
 388          - lat_1: Latitude of first standard parallel
 389          - lat_2: Latitude of second standard parallel
 390         \par Example
 391         \image html ex_euler.gif
 392     */
 393     template <typename T, typename Parameters>
 394     struct euler_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
 395     {
 396         template <typename Params>
 397         inline euler_spheroid(Params const& params, Parameters& par)
 398         {
 399             detail::sconics::setup_euler(params, par, this->m_proj_parm);
 400         }
 401     };
 402
 403     /*!
 404         \brief Perspective Conic projection
 405         \ingroup projections
 406         \tparam Geographic latlong point type
 407         \tparam Cartesian xy point type
 408         \tparam Parameters parameter type
 409         \par Projection characteristics
 410          - Conic
 411          - Spheroid
 412         \par Projection parameters
 413          - lat_1: Latitude of first standard parallel
 414          - lat_2: Latitude of second standard parallel
 415         \par Example
 416         \image html ex_pconic.gif
 417     */
 418     template <typename T, typename Parameters>
 419     struct pconic_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
 420     {
 421         template <typename Params>
 422         inline pconic_spheroid(Params const& params, Parameters& par)
 423         {
 424             detail::sconics::setup_pconic(params, par, this->m_proj_parm);
 425         }
 426     };
 427
 428     /*!
 429         \brief Vitkovsky I projection
 430         \ingroup projections
 431         \tparam Geographic latlong point type
 432         \tparam Cartesian xy point type
 433         \tparam Parameters parameter type
 434         \par Projection characteristics
 435          - Conic
 436          - Spheroid
 437         \par Projection parameters
 438          - lat_1: Latitude of first standard parallel
 439          - lat_2: Latitude of second standard parallel
 440         \par Example
 441         \image html ex_vitk1.gif
 442     */
 443     template <typename T, typename Parameters>
 444     struct vitk1_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
 445     {
 446         template <typename Params>
 447         inline vitk1_spheroid(Params const& params, Parameters& par)
 448         {
 449             detail::sconics::setup_vitk1(params, par, this->m_proj_parm);
 450         }
 451     };
 452
 453     #ifndef DOXYGEN_NO_DETAIL
 454     namespace detail
 455     {
 456
 457         // Static projection
 458         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_euler, euler_spheroid)
 459         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_murd1, murd1_spheroid)
 460         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_murd2, murd2_spheroid)
 461         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_murd3, murd3_spheroid)
 462         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_pconic, pconic_spheroid)
 463         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_tissot, tissot_spheroid)
 464         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_vitk1, vitk1_spheroid)
 465
 466         // Factory entry(s)
 467         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(euler_entry, euler_spheroid)
 468         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(murd1_entry, murd1_spheroid)
 469         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(murd2_entry, murd2_spheroid)
 470         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(murd3_entry, murd3_spheroid)
 471         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(pconic_entry, pconic_spheroid)
 472         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(tissot_entry, tissot_spheroid)
 473         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(vitk1_entry, vitk1_spheroid)
 474
 475         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(sconics_init)
 476         {
 477             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(euler, euler_entry)
 478             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(murd1, murd1_entry)
 479             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(murd2, murd2_entry)
 480             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(murd3, murd3_entry)
 481             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(pconic, pconic_entry)
 482             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(tissot, tissot_entry)
 483             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(vitk1, vitk1_entry)
 484         }
 485
 486     } // namespace detail
 487     #endif // doxygen
 488
 489 } // namespace projections
 490
 491 }} // namespace boost::geometry
 492
 493 #endif // BOOST_GEOMETRY_PROJECTIONS_SCONICS_HPP
 494