boost/geometry/srs/projections/proj/mod_ster.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_MOD_STER_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
  42
  43 #include <boost/geometry/util/math.hpp>
  44 #include <boost/math/special_functions/hypot.hpp>
  45
  46 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  47 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  48 #include <boost/geometry/srs/projections/impl/projects.hpp>
  49 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  50 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
  51 #include <boost/geometry/srs/projections/impl/pj_zpoly1.hpp>
  52
  53 namespace boost { namespace geometry
  54 {
  55
  56 namespace projections
  57 {
  58     #ifndef DOXYGEN_NO_DETAIL
  59     namespace detail { namespace mod_ster
  60     {
  61
  62             static const double epsilon = 1e-12;
  63
  64             template <typename T>
  65             struct par_mod_ster
  66             {
  67                 T              cchio, schio;
  68                 pj_complex<T>* zcoeff;
  69                 int            n;
  70             };
  71
  72             /* based upon Snyder and Linck, USGS-NMD */
  73
  74             template <typename T, typename Parameters>
  75             struct base_mod_ster_ellipsoid
  76             {
  77                 par_mod_ster<T> m_proj_parm;
  78
  79                 // FORWARD(e_forward)  ellipsoid
  80                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  81                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  82                 {
  83                     static const T half_pi = detail::half_pi<T>();
  84
  85                     T sinlon, coslon, esphi, chi, schi, cchi, s;
  86                     pj_complex<T> p;
  87
  88                     sinlon = sin(lp_lon);
  89                     coslon = cos(lp_lon);
  90                     esphi = par.e * sin(lp_lat);
  91                     chi = 2. * atan(tan((half_pi + lp_lat) * .5) *
  92                         math::pow((T(1) - esphi) / (T(1) + esphi), par.e * T(0.5))) - half_pi;
  93                     schi = sin(chi);
  94                     cchi = cos(chi);
  95                     s = 2. / (1. + this->m_proj_parm.schio * schi + this->m_proj_parm.cchio * cchi * coslon);
  96                     p.r = s * cchi * sinlon;
  97                     p.i = s * (this->m_proj_parm.cchio * schi - this->m_proj_parm.schio * cchi * coslon);
  98                     p = pj_zpoly1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n);
  99                     xy_x = p.r;
 100                     xy_y = p.i;
 101                 }
 102
 103                 // INVERSE(e_inverse)  ellipsoid
 104                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 105                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
 106                 {
 107                     static const T half_pi = detail::half_pi<T>();
 108
 109                     int nn;
 110                     pj_complex<T> p, fxy, fpxy, dp;
 111                     T den, rh = 0, z, sinz = 0, cosz = 0, chi, phi = 0, dphi, esphi;
 112
 113                     p.r = xy_x;
 114                     p.i = xy_y;
 115                     for (nn = 20; nn ;--nn) {
 116                         fxy = pj_zpolyd1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n, &fpxy);
 117                         fxy.r -= xy_x;
 118                         fxy.i -= xy_y;
 119                         den = fpxy.r * fpxy.r + fpxy.i * fpxy.i;
 120                         dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den;
 121                         dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den;
 122                         p.r += dp.r;
 123                         p.i += dp.i;
 124                         if ((fabs(dp.r) + fabs(dp.i)) <= epsilon)
 125                             break;
 126                     }
 127                     if (nn) {
 128                         rh = boost::math::hypot(p.r, p.i);
 129                         z = 2. * atan(.5 * rh);
 130                         sinz = sin(z);
 131                         cosz = cos(z);
 132                         lp_lon = par.lam0;
 133                         if (fabs(rh) <= epsilon) {
 134                             /* if we end up here input coordinates were (0,0).
 135                              * pj_inv() adds P->lam0 to lp.lam, this way we are
 136                              * sure to get the correct offset */
 137                             lp_lon = 0.0;
 138                             lp_lat = par.phi0;
 139                             return;
 140                         }
 141                         chi = aasin(cosz * this->m_proj_parm.schio + p.i * sinz * this->m_proj_parm.cchio / rh);
 142                         phi = chi;
 143                         for (nn = 20; nn ;--nn) {
 144                             esphi = par.e * sin(phi);
 145                             dphi = 2. * atan(tan((half_pi + chi) * .5) *
 146                                 math::pow((T(1) + esphi) / (T(1) - esphi), par.e * T(0.5))) - half_pi - phi;
 147                             phi += dphi;
 148                             if (fabs(dphi) <= epsilon)
 149                                 break;
 150                         }
 151                     }
 152                     if (nn) {
 153                         lp_lat = phi;
 154                         lp_lon = atan2(p.r * sinz, rh * this->m_proj_parm.cchio * cosz - p.i *
 155                             this->m_proj_parm.schio * sinz);
 156                     } else
 157                         lp_lon = lp_lat = HUGE_VAL;
 158                 }
 159
 160                 static inline std::string get_name()
 161                 {
 162                     return "mod_ster_ellipsoid";
 163                 }
 164
 165             };
 166
 167             template <typename Parameters, typename T>
 168             inline void setup(Parameters& par, par_mod_ster<T>& proj_parm)  /* general initialization */
 169             {
 170                 static T const half_pi = detail::half_pi<T>();
 171
 172                 T esphi, chio;
 173
 174                 if (par.es != 0.0) {
 175                     esphi = par.e * sin(par.phi0);
 176                     chio = 2. * atan(tan((half_pi + par.phi0) * .5) *
 177                         math::pow((T(1) - esphi) / (T(1) + esphi), par.e * T(0.5))) - half_pi;
 178                 } else
 179                     chio = par.phi0;
 180                 proj_parm.schio = sin(chio);
 181                 proj_parm.cchio = cos(chio);
 182             }
 183
 184
 185             /* Miller Oblated Stereographic */
 186             template <typename Parameters, typename T>
 187             inline void setup_mil_os(Parameters& par, par_mod_ster<T>& proj_parm)
 188             {
 189                 static const T d2r = geometry::math::d2r<T>();
 190
 191                 static pj_complex<T> AB[] = {
 192                     {0.924500, 0.},
 193                     {0.,       0.},
 194                     {0.019430, 0.}
 195                 };
 196
 197                 proj_parm.n = 2;
 198                 par.lam0 = d2r * 20.;
 199                 par.phi0 = d2r * 18.;
 200                 proj_parm.zcoeff = AB;
 201                 par.es = 0.;
 202
 203                 setup(par, proj_parm);
 204             }
 205
 206             /* Lee Oblated Stereographic */
 207             template <typename Parameters, typename T>
 208             inline void setup_lee_os(Parameters& par, par_mod_ster<T>& proj_parm)
 209             {
 210                 static const T d2r = geometry::math::d2r<T>();
 211
 212                 static pj_complex<T> AB[] = {
 213                     { 0.721316,   0.},
 214                     { 0.,         0.},
 215                     {-0.0088162, -0.00617325}
 216                 };
 217
 218                 proj_parm.n = 2;
 219                 par.lam0 = d2r * -165.;
 220                 par.phi0 = d2r * -10.;
 221                 proj_parm.zcoeff = AB;
 222                 par.es = 0.;
 223
 224                 setup(par, proj_parm);
 225             }
 226
 227             // Mod. Stererographics of 48 U.S.
 228             template <typename Parameters, typename T>
 229             inline void setup_gs48(Parameters& par, par_mod_ster<T>& proj_parm)
 230             {
 231                 static const T d2r = geometry::math::d2r<T>();
 232
 233                 static pj_complex<T> AB[] = { /* 48 United States */
 234                     { 0.98879,  0.},
 235                     { 0.,       0.},
 236                     {-0.050909, 0.},
 237                     { 0.,       0.},
 238                     { 0.075528, 0.}
 239                 };
 240
 241                 proj_parm.n = 4;
 242                 par.lam0 = d2r * -96.;
 243                 par.phi0 = d2r * -39.;
 244                 proj_parm.zcoeff = AB;
 245                 par.es = 0.;
 246                 par.a = 6370997.;
 247
 248                 setup(par, proj_parm);
 249             }
 250
 251             // Mod. Stererographics of Alaska
 252             template <typename Parameters, typename T>
 253             inline void setup_alsk(Parameters& par, par_mod_ster<T>& proj_parm)
 254             {
 255                 static const T d2r = geometry::math::d2r<T>();
 256
 257                 static pj_complex<T> ABe[] = { /* Alaska ellipsoid */
 258                     { .9945303, 0.},
 259                     { .0052083, -.0027404},
 260                     { .0072721,  .0048181},
 261                     {-.0151089, -.1932526},
 262                     { .0642675, -.1381226},
 263                     { .3582802, -.2884586}
 264                 };
 265
 266                 static pj_complex<T> ABs[] = { /* Alaska sphere */
 267                     { .9972523, 0.},
 268                     { .0052513, -.0041175},
 269                     { .0074606,  .0048125},
 270                     {-.0153783, -.1968253},
 271                     { .0636871, -.1408027},
 272                     { .3660976, -.2937382}
 273                 };
 274
 275                 proj_parm.n = 5;
 276                 par.lam0 = d2r * -152.;
 277                 par.phi0 = d2r * 64.;
 278                 if (par.es != 0.0) { /* fixed ellipsoid/sphere */
 279                     proj_parm.zcoeff = ABe;
 280                     par.a = 6378206.4;
 281                     par.e = sqrt(par.es = 0.00676866);
 282                 } else {
 283                     proj_parm.zcoeff = ABs;
 284                     par.a = 6370997.;
 285                 }
 286
 287                 setup(par, proj_parm);
 288             }
 289
 290             // Mod. Stererographics of 50 U.S.
 291             template <typename Parameters, typename T>
 292             inline void setup_gs50(Parameters& par, par_mod_ster<T>& proj_parm)
 293             {
 294                 static const T d2r = geometry::math::d2r<T>();
 295
 296                 static pj_complex<T> ABe[] = { /* GS50 ellipsoid */
 297                     { .9827497, 0.},
 298                     { .0210669,  .0053804},
 299                     {-.1031415, -.0571664},
 300                     {-.0323337, -.0322847},
 301                     { .0502303,  .1211983},
 302                     { .0251805,  .0895678},
 303                     {-.0012315, -.1416121},
 304                     { .0072202, -.1317091},
 305                     {-.0194029,  .0759677},
 306                     {-.0210072,  .0834037}
 307                 };
 308                 static pj_complex<T> ABs[] = { /* GS50 sphere */
 309                     { .9842990, 0.},
 310                     { .0211642,  .0037608},
 311                     {-.1036018, -.0575102},
 312                     {-.0329095, -.0320119},
 313                     { .0499471,  .1223335},
 314                     { .0260460,  .0899805},
 315                     { .0007388, -.1435792},
 316                     { .0075848, -.1334108},
 317                     {-.0216473,  .0776645},
 318                     {-.0225161,  .0853673}
 319                 };
 320
 321                 proj_parm.n = 9;
 322                 par.lam0 = d2r * -120.;
 323                 par.phi0 = d2r * 45.;
 324                 if (par.es != 0.0) { /* fixed ellipsoid/sphere */
 325                     proj_parm.zcoeff = ABe;
 326                     par.a = 6378206.4;
 327                     par.e = sqrt(par.es = 0.00676866);
 328                 } else {
 329                     proj_parm.zcoeff = ABs;
 330                     par.a = 6370997.;
 331                 }
 332
 333                 setup(par, proj_parm);
 334             }
 335
 336     }} // namespace detail::mod_ster
 337     #endif // doxygen
 338
 339     /*!
 340         \brief Miller Oblated Stereographic projection
 341         \ingroup projections
 342         \tparam Geographic latlong point type
 343         \tparam Cartesian xy point type
 344         \tparam Parameters parameter type
 345         \par Projection characteristics
 346          - Azimuthal (mod)
 347         \par Example
 348         \image html ex_mil_os.gif
 349     */
 350     template <typename T, typename Parameters>
 351     struct mil_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
 352     {
 353         template <typename Params>
 354         inline mil_os_ellipsoid(Params const& , Parameters & par)
 355         {
 356             detail::mod_ster::setup_mil_os(par, this->m_proj_parm);
 357         }
 358     };
 359
 360     /*!
 361         \brief Lee Oblated Stereographic projection
 362         \ingroup projections
 363         \tparam Geographic latlong point type
 364         \tparam Cartesian xy point type
 365         \tparam Parameters parameter type
 366         \par Projection characteristics
 367          - Azimuthal (mod)
 368         \par Example
 369         \image html ex_lee_os.gif
 370     */
 371     template <typename T, typename Parameters>
 372     struct lee_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
 373     {
 374         template <typename Params>
 375         inline lee_os_ellipsoid(Params const& , Parameters & par)
 376         {
 377             detail::mod_ster::setup_lee_os(par, this->m_proj_parm);
 378         }
 379     };
 380
 381     /*!
 382         \brief Mod. Stererographics of 48 U.S. projection
 383         \ingroup projections
 384         \tparam Geographic latlong point type
 385         \tparam Cartesian xy point type
 386         \tparam Parameters parameter type
 387         \par Projection characteristics
 388          - Azimuthal (mod)
 389         \par Example
 390         \image html ex_gs48.gif
 391     */
 392     template <typename T, typename Parameters>
 393     struct gs48_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
 394     {
 395         template <typename Params>
 396         inline gs48_ellipsoid(Params const& , Parameters & par)
 397         {
 398             detail::mod_ster::setup_gs48(par, this->m_proj_parm);
 399         }
 400     };
 401
 402     /*!
 403         \brief Mod. Stererographics of Alaska projection
 404         \ingroup projections
 405         \tparam Geographic latlong point type
 406         \tparam Cartesian xy point type
 407         \tparam Parameters parameter type
 408         \par Projection characteristics
 409          - Azimuthal (mod)
 410         \par Example
 411         \image html ex_alsk.gif
 412     */
 413     template <typename T, typename Parameters>
 414     struct alsk_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
 415     {
 416         template <typename Params>
 417         inline alsk_ellipsoid(Params const& , Parameters & par)
 418         {
 419             detail::mod_ster::setup_alsk(par, this->m_proj_parm);
 420         }
 421     };
 422
 423     /*!
 424         \brief Mod. Stererographics of 50 U.S. projection
 425         \ingroup projections
 426         \tparam Geographic latlong point type
 427         \tparam Cartesian xy point type
 428         \tparam Parameters parameter type
 429         \par Projection characteristics
 430          - Azimuthal (mod)
 431         \par Example
 432         \image html ex_gs50.gif
 433     */
 434     template <typename T, typename Parameters>
 435     struct gs50_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
 436     {
 437         template <typename Params>
 438         inline gs50_ellipsoid(Params const& , Parameters & par)
 439         {
 440             detail::mod_ster::setup_gs50(par, this->m_proj_parm);
 441         }
 442     };
 443
 444     #ifndef DOXYGEN_NO_DETAIL
 445     namespace detail
 446     {
 447
 448         // Static projection
 449         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_mil_os, mil_os_ellipsoid)
 450         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lee_os, lee_os_ellipsoid)
 451         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_gs48, gs48_ellipsoid)
 452         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_alsk, alsk_ellipsoid)
 453         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_gs50, gs50_ellipsoid)
 454
 455         // Factory entry(s)
 456         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(mil_os_entry, mil_os_ellipsoid)
 457         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lee_os_entry, lee_os_ellipsoid)
 458         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(gs48_entry, gs48_ellipsoid)
 459         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(alsk_entry, alsk_ellipsoid)
 460         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(gs50_entry, gs50_ellipsoid)
 461
 462         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(mod_ster_init)
 463         {
 464             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(mil_os, mil_os_entry)
 465             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lee_os, lee_os_entry)
 466             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(gs48, gs48_entry)
 467             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(alsk, alsk_entry)
 468             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(gs50, gs50_entry)
 469         }
 470
 471     } // namespace detail
 472     #endif // doxygen
 473
 474 } // namespace projections
 475
 476 }} // namespace boost::geometry
 477
 478 #endif // BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
 479