boost/geometry/srs/projections/proj/stere.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_STERE_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
  42
  43 #include <boost/config.hpp>
  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/pj_tsfn.hpp>
  52 #include <boost/geometry/srs/projections/impl/projects.hpp>
  53
  54 namespace boost { namespace geometry
  55 {
  56
  57 namespace projections
  58 {
  59     #ifndef DOXYGEN_NO_DETAIL
  60     namespace detail { namespace stere
  61     {
  62             static const double epsilon10 = 1.e-10;
  63             static const double tolerance = 1.e-8;
  64             static const int n_iter = 8;
  65             static const double conv_tolerance = 1.e-10;
  66
  67             enum mode_type {
  68                 s_pole = 0,
  69                 n_pole = 1,
  70                 obliq  = 2,
  71                 equit  = 3
  72             };
  73
  74             template <typename T>
  75             struct par_stere
  76             {
  77                 T   phits;
  78                 T   sinX1;
  79                 T   cosX1;
  80                 T   akm1;
  81                 mode_type mode;
  82             };
  83
  84             template <typename T>
  85             inline T ssfn_(T const& phit, T sinphi, T const& eccen)
  86             {
  87                 static const T half_pi = detail::half_pi<T>();
  88
  89                 sinphi *= eccen;
  90                 return (tan (.5 * (half_pi + phit)) *
  91                    math::pow((T(1) - sinphi) / (T(1) + sinphi), T(0.5) * eccen));
  92             }
  93
  94             template <typename T, typename Parameters>
  95             struct base_stere_ellipsoid
  96             {
  97                 par_stere<T> m_proj_parm;
  98
  99                 // FORWARD(e_forward)  ellipsoid
 100                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 101                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
 102                 {
 103                     static const T half_pi = detail::half_pi<T>();
 104
 105                     T coslam, sinlam, sinX=0.0, cosX=0.0, X, A = 0.0, sinphi;
 106
 107                     coslam = cos(lp_lon);
 108                     sinlam = sin(lp_lon);
 109                     sinphi = sin(lp_lat);
 110                     if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
 111                         sinX = sin(X = 2. * atan(ssfn_(lp_lat, sinphi, par.e)) - half_pi);
 112                         cosX = cos(X);
 113                     }
 114                     switch (this->m_proj_parm.mode) {
 115                     case obliq:
 116                         A = this->m_proj_parm.akm1 / (this->m_proj_parm.cosX1 * (1. + this->m_proj_parm.sinX1 * sinX +
 117                            this->m_proj_parm.cosX1 * cosX * coslam));
 118                         xy_y = A * (this->m_proj_parm.cosX1 * sinX - this->m_proj_parm.sinX1 * cosX * coslam);
 119                         goto xmul; /* but why not just  xy.x = A * cosX; break;  ? */
 120
 121                     case equit:
 122                         // TODO: calculate denominator once
 123                         /* avoid zero division */
 124                         if (1. + cosX * coslam == 0.0) {
 125                             xy_y = HUGE_VAL;
 126                         } else {
 127                             A = this->m_proj_parm.akm1 / (1. + cosX * coslam);
 128                             xy_y = A * sinX;
 129                         }
 130                 xmul:
 131                         xy_x = A * cosX;
 132                         break;
 133
 134                     case s_pole:
 135                         lp_lat = -lp_lat;
 136                         coslam = - coslam;
 137                         sinphi = -sinphi;
 138                         BOOST_FALLTHROUGH;
 139                     case n_pole:
 140                         xy_x = this->m_proj_parm.akm1 * pj_tsfn(lp_lat, sinphi, par.e);
 141                         xy_y = - xy_x * coslam;
 142                         break;
 143                     }
 144
 145                     xy_x = xy_x * sinlam;
 146                 }
 147
 148                 // INVERSE(e_inverse)  ellipsoid
 149                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 150                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 151                 {
 152                     static const T half_pi = detail::half_pi<T>();
 153
 154                     T cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, halfpi=0.0;
 155                     int i;
 156
 157                     rho = boost::math::hypot(xy_x, xy_y);
 158                     switch (this->m_proj_parm.mode) {
 159                     case obliq:
 160                     case equit:
 161                         cosphi = cos( tp = 2. * atan2(rho * this->m_proj_parm.cosX1 , this->m_proj_parm.akm1) );
 162                         sinphi = sin(tp);
 163                         if( rho == 0.0 )
 164                             phi_l = asin(cosphi * this->m_proj_parm.sinX1);
 165                         else
 166                             phi_l = asin(cosphi * this->m_proj_parm.sinX1 + (xy_y * sinphi * this->m_proj_parm.cosX1 / rho));
 167
 168                         tp = tan(.5 * (half_pi + phi_l));
 169                         xy_x *= sinphi;
 170                         xy_y = rho * this->m_proj_parm.cosX1 * cosphi - xy_y * this->m_proj_parm.sinX1* sinphi;
 171                         halfpi = half_pi;
 172                         halfe = .5 * par.e;
 173                         break;
 174                     case n_pole:
 175                         xy_y = -xy_y;
 176                         BOOST_FALLTHROUGH;
 177                     case s_pole:
 178                         phi_l = half_pi - 2. * atan(tp = - rho / this->m_proj_parm.akm1);
 179                         halfpi = -half_pi;
 180                         halfe = -.5 * par.e;
 181                         break;
 182                     }
 183                     for (i = n_iter; i--; phi_l = lp_lat) {
 184                         sinphi = par.e * sin(phi_l);
 185                         lp_lat = T(2) * atan(tp * math::pow((T(1)+sinphi)/(T(1)-sinphi), halfe)) - halfpi;
 186                         if (fabs(phi_l - lp_lat) < conv_tolerance) {
 187                             if (this->m_proj_parm.mode == s_pole)
 188                                 lp_lat = -lp_lat;
 189                             lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
 190                             return;
 191                         }
 192                     }
 193                     BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 194                 }
 195
 196                 static inline std::string get_name()
 197                 {
 198                     return "stere_ellipsoid";
 199                 }
 200
 201             };
 202
 203             template <typename T, typename Parameters>
 204             struct base_stere_spheroid
 205             {
 206                 par_stere<T> m_proj_parm;
 207
 208                 // FORWARD(s_forward)  spheroid
 209                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 210                 inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
 211                 {
 212                     static const T fourth_pi = detail::fourth_pi<T>();
 213                     static const T half_pi = detail::half_pi<T>();
 214
 215                     T  sinphi, cosphi, coslam, sinlam;
 216
 217                     sinphi = sin(lp_lat);
 218                     cosphi = cos(lp_lat);
 219                     coslam = cos(lp_lon);
 220                     sinlam = sin(lp_lon);
 221                     switch (this->m_proj_parm.mode) {
 222                     case equit:
 223                         xy_y = 1. + cosphi * coslam;
 224                         goto oblcon;
 225                     case obliq:
 226                         xy_y = 1. + this->m_proj_parm.sinX1 * sinphi + this->m_proj_parm.cosX1 * cosphi * coslam;
 227                 oblcon:
 228                         if (xy_y <= epsilon10) {
 229                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 230                         }
 231                         xy_x = (xy_y = this->m_proj_parm.akm1 / xy_y) * cosphi * sinlam;
 232                         xy_y *= (this->m_proj_parm.mode == equit) ? sinphi :
 233                            this->m_proj_parm.cosX1 * sinphi - this->m_proj_parm.sinX1 * cosphi * coslam;
 234                         break;
 235                     case n_pole:
 236                         coslam = - coslam;
 237                         lp_lat = - lp_lat;
 238                         BOOST_FALLTHROUGH;
 239                     case s_pole:
 240                         if (fabs(lp_lat - half_pi) < tolerance) {
 241                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 242                         }
 243                         xy_x = sinlam * ( xy_y = this->m_proj_parm.akm1 * tan(fourth_pi + .5 * lp_lat) );
 244                         xy_y *= coslam;
 245                         break;
 246                     }
 247                 }
 248
 249                 // INVERSE(s_inverse)  spheroid
 250                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 251                 inline void inv(Parameters const& par, T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 252                 {
 253                     T  c, rh, sinc, cosc;
 254
 255                     sinc = sin(c = 2. * atan((rh = boost::math::hypot(xy_x, xy_y)) / this->m_proj_parm.akm1));
 256                     cosc = cos(c);
 257                     lp_lon = 0.;
 258
 259                     switch (this->m_proj_parm.mode) {
 260                     case equit:
 261                         if (fabs(rh) <= epsilon10)
 262                             lp_lat = 0.;
 263                         else
 264                             lp_lat = asin(xy_y * sinc / rh);
 265                         if (cosc != 0. || xy_x != 0.)
 266                             lp_lon = atan2(xy_x * sinc, cosc * rh);
 267                         break;
 268                     case obliq:
 269                         if (fabs(rh) <= epsilon10)
 270                             lp_lat = par.phi0;
 271                         else
 272                             lp_lat = asin(cosc * this->m_proj_parm.sinX1 + xy_y * sinc * this->m_proj_parm.cosX1 / rh);
 273                         if ((c = cosc - this->m_proj_parm.sinX1 * sin(lp_lat)) != 0. || xy_x != 0.)
 274                             lp_lon = atan2(xy_x * sinc * this->m_proj_parm.cosX1, c * rh);
 275                         break;
 276                     case n_pole:
 277                         xy_y = -xy_y;
 278                         BOOST_FALLTHROUGH;
 279                     case s_pole:
 280                         if (fabs(rh) <= epsilon10)
 281                             lp_lat = par.phi0;
 282                         else
 283                             lp_lat = asin(this->m_proj_parm.mode == s_pole ? - cosc : cosc);
 284                         lp_lon = (xy_x == 0. && xy_y == 0.) ? 0. : atan2(xy_x, xy_y);
 285                         break;
 286                     }
 287                 }
 288
 289                 static inline std::string get_name()
 290                 {
 291                     return "stere_spheroid";
 292                 }
 293
 294             };
 295
 296             template <typename Parameters, typename T>
 297             inline void setup(Parameters const& par, par_stere<T>& proj_parm)  /* general initialization */
 298             {
 299                 static const T fourth_pi = detail::fourth_pi<T>();
 300                 static const T half_pi = detail::half_pi<T>();
 301
 302                 T t;
 303
 304                 if (fabs((t = fabs(par.phi0)) - half_pi) < epsilon10)
 305                     proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
 306                 else
 307                     proj_parm.mode = t > epsilon10 ? obliq : equit;
 308                 proj_parm.phits = fabs(proj_parm.phits);
 309
 310                 if (par.es != 0.0) {
 311                     T X;
 312
 313                     switch (proj_parm.mode) {
 314                     case n_pole:
 315                     case s_pole:
 316                         if (fabs(proj_parm.phits - half_pi) < epsilon10)
 317                             proj_parm.akm1 = 2. * par.k0 /
 318                                sqrt(math::pow(T(1)+par.e,T(1)+par.e)*math::pow(T(1)-par.e,T(1)-par.e));
 319                         else {
 320                             proj_parm.akm1 = cos(proj_parm.phits) /
 321                                pj_tsfn(proj_parm.phits, t = sin(proj_parm.phits), par.e);
 322                             t *= par.e;
 323                             proj_parm.akm1 /= sqrt(1. - t * t);
 324                         }
 325                         break;
 326                     case equit:
 327                     case obliq:
 328                         t = sin(par.phi0);
 329                         X = 2. * atan(ssfn_(par.phi0, t, par.e)) - half_pi;
 330                         t *= par.e;
 331                         proj_parm.akm1 = 2. * par.k0 * cos(par.phi0) / sqrt(1. - t * t);
 332                         proj_parm.sinX1 = sin(X);
 333                         proj_parm.cosX1 = cos(X);
 334                         break;
 335                     }
 336                 } else {
 337                     switch (proj_parm.mode) {
 338                     case obliq:
 339                         proj_parm.sinX1 = sin(par.phi0);
 340                         proj_parm.cosX1 = cos(par.phi0);
 341                         BOOST_FALLTHROUGH;
 342                     case equit:
 343                         proj_parm.akm1 = 2. * par.k0;
 344                         break;
 345                     case s_pole:
 346                     case n_pole:
 347                         proj_parm.akm1 = fabs(proj_parm.phits - half_pi) >= epsilon10 ?
 348                            cos(proj_parm.phits) / tan(fourth_pi - .5 * proj_parm.phits) :
 349                            2. * par.k0 ;
 350                         break;
 351                     }
 352                 }
 353             }
 354
 355
 356             // Stereographic
 357             template <typename Params, typename Parameters, typename T>
 358             inline void setup_stere(Params const& params, Parameters const& par, par_stere<T>& proj_parm)
 359             {
 360                 static const T half_pi = detail::half_pi<T>();
 361
 362                 if (! pj_param_r<srs::spar::lat_ts>(params, "lat_ts", srs::dpar::lat_ts, proj_parm.phits))
 363                     proj_parm.phits = half_pi;
 364
 365                 setup(par, proj_parm);
 366             }
 367
 368             // Universal Polar Stereographic
 369             template <typename Params, typename Parameters, typename T>
 370             inline void setup_ups(Params const& params, Parameters& par, par_stere<T>& proj_parm)
 371             {
 372                 static const T half_pi = detail::half_pi<T>();
 373
 374                 /* International Ellipsoid */
 375                 par.phi0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? -half_pi: half_pi;
 376                 if (par.es == 0.0) {
 377                     BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
 378                 }
 379                 par.k0 = .994;
 380                 par.x0 = 2000000.;
 381                 par.y0 = 2000000.;
 382                 proj_parm.phits = half_pi;
 383                 par.lam0 = 0.;
 384
 385                 setup(par, proj_parm);
 386             }
 387
 388     }} // namespace detail::stere
 389     #endif // doxygen
 390
 391     /*!
 392         \brief Stereographic projection
 393         \ingroup projections
 394         \tparam Geographic latlong point type
 395         \tparam Cartesian xy point type
 396         \tparam Parameters parameter type
 397         \par Projection characteristics
 398          - Azimuthal
 399          - Spheroid
 400          - Ellipsoid
 401         \par Projection parameters
 402          - lat_ts: Latitude of true scale (degrees)
 403         \par Example
 404         \image html ex_stere.gif
 405     */
 406     template <typename T, typename Parameters>
 407     struct stere_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
 408     {
 409         template <typename Params>
 410         inline stere_ellipsoid(Params const& params, Parameters const& par)
 411         {
 412             detail::stere::setup_stere(params, par, this->m_proj_parm);
 413         }
 414     };
 415
 416     /*!
 417         \brief Stereographic projection
 418         \ingroup projections
 419         \tparam Geographic latlong point type
 420         \tparam Cartesian xy point type
 421         \tparam Parameters parameter type
 422         \par Projection characteristics
 423          - Azimuthal
 424          - Spheroid
 425          - Ellipsoid
 426         \par Projection parameters
 427          - lat_ts: Latitude of true scale (degrees)
 428         \par Example
 429         \image html ex_stere.gif
 430     */
 431     template <typename T, typename Parameters>
 432     struct stere_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
 433     {
 434         template <typename Params>
 435         inline stere_spheroid(Params const& params, Parameters const& par)
 436         {
 437             detail::stere::setup_stere(params, par, this->m_proj_parm);
 438         }
 439     };
 440
 441     /*!
 442         \brief Universal Polar Stereographic projection
 443         \ingroup projections
 444         \tparam Geographic latlong point type
 445         \tparam Cartesian xy point type
 446         \tparam Parameters parameter type
 447         \par Projection characteristics
 448          - Azimuthal
 449          - Spheroid
 450          - Ellipsoid
 451         \par Projection parameters
 452          - south: Denotes southern hemisphere UTM zone (boolean)
 453         \par Example
 454         \image html ex_ups.gif
 455     */
 456     template <typename T, typename Parameters>
 457     struct ups_ellipsoid : public detail::stere::base_stere_ellipsoid<T, Parameters>
 458     {
 459         template <typename Params>
 460         inline ups_ellipsoid(Params const& params, Parameters & par)
 461         {
 462             detail::stere::setup_ups(params, par, this->m_proj_parm);
 463         }
 464     };
 465
 466     /*!
 467         \brief Universal Polar Stereographic projection
 468         \ingroup projections
 469         \tparam Geographic latlong point type
 470         \tparam Cartesian xy point type
 471         \tparam Parameters parameter type
 472         \par Projection characteristics
 473          - Azimuthal
 474          - Spheroid
 475          - Ellipsoid
 476         \par Projection parameters
 477          - south: Denotes southern hemisphere UTM zone (boolean)
 478         \par Example
 479         \image html ex_ups.gif
 480     */
 481     template <typename T, typename Parameters>
 482     struct ups_spheroid : public detail::stere::base_stere_spheroid<T, Parameters>
 483     {
 484         template <typename Params>
 485         inline ups_spheroid(Params const& params, Parameters & par)
 486         {
 487             detail::stere::setup_ups(params, par, this->m_proj_parm);
 488         }
 489     };
 490
 491     #ifndef DOXYGEN_NO_DETAIL
 492     namespace detail
 493     {
 494
 495         // Static projection
 496         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_stere, stere_spheroid, stere_ellipsoid)
 497         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_ups, ups_spheroid, ups_ellipsoid)
 498
 499         // Factory entry(s)
 500         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(stere_entry, stere_spheroid, stere_ellipsoid)
 501         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(ups_entry, ups_spheroid, ups_ellipsoid)
 502
 503         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(stere_init)
 504         {
 505             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(stere, stere_entry)
 506             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(ups, ups_entry)
 507         }
 508
 509     } // namespace detail
 510     #endif // doxygen
 511
 512 } // namespace projections
 513
 514 }} // namespace boost::geometry
 515
 516 #endif // BOOST_GEOMETRY_PROJECTIONS_STERE_HPP
 517