boost/geometry/srs/projections/proj/laea.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_LAEA_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_LAEA_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/projects.hpp>
  50 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  51 #include <boost/geometry/srs/projections/impl/pj_auth.hpp>
  52 #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
  53
  54 namespace boost { namespace geometry
  55 {
  56
  57 namespace projections
  58 {
  59     #ifndef DOXYGEN_NO_DETAIL
  60     namespace detail { namespace laea
  61     {
  62             static const double epsilon10 = 1.e-10;
  63
  64             enum mode_type {
  65                 n_pole = 0,
  66                 s_pole = 1,
  67                 equit  = 2,
  68                 obliq  = 3
  69             };
  70
  71             template <typename T>
  72             struct par_laea
  73             {
  74                 T   sinb1;
  75                 T   cosb1;
  76                 T   xmf;
  77                 T   ymf;
  78                 T   mmf;
  79                 T   qp;
  80                 T   dd;
  81                 T   rq;
  82                 detail::apa<T> apa;
  83                 mode_type mode;
  84             };
  85
  86             template <typename T, typename Parameters>
  87             struct base_laea_ellipsoid
  88             {
  89                 par_laea<T> m_proj_parm;
  90
  91                 // FORWARD(e_forward)  ellipsoid
  92                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  93                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  94                 {
  95                     static const T half_pi = detail::half_pi<T>();
  96
  97                     T coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
  98
  99                     coslam = cos(lp_lon);
 100                     sinlam = sin(lp_lon);
 101                     sinphi = sin(lp_lat);
 102                     q = pj_qsfn(sinphi, par.e, par.one_es);
 103
 104                     if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
 105                         sinb = q / this->m_proj_parm.qp;
 106                         cosb = sqrt(1. - sinb * sinb);
 107                     }
 108
 109                     switch (this->m_proj_parm.mode) {
 110                     case obliq:
 111                         b = 1. + this->m_proj_parm.sinb1 * sinb + this->m_proj_parm.cosb1 * cosb * coslam;
 112                         break;
 113                     case equit:
 114                         b = 1. + cosb * coslam;
 115                         break;
 116                     case n_pole:
 117                         b = half_pi + lp_lat;
 118                         q = this->m_proj_parm.qp - q;
 119                         break;
 120                     case s_pole:
 121                         b = lp_lat - half_pi;
 122                         q = this->m_proj_parm.qp + q;
 123                         break;
 124                     }
 125                     if (fabs(b) < epsilon10) {
 126                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 127                     }
 128
 129                     switch (this->m_proj_parm.mode) {
 130                     case obliq:
 131                         b = sqrt(2. / b);
 132                         xy_y = this->m_proj_parm.ymf * b * (this->m_proj_parm.cosb1 * sinb - this->m_proj_parm.sinb1 * cosb * coslam);
 133                         goto eqcon;
 134                         break;
 135                     case equit:
 136                         b = sqrt(2. / (1. + cosb * coslam));
 137                         xy_y = b * sinb * this->m_proj_parm.ymf;
 138                 eqcon:
 139                         xy_x = this->m_proj_parm.xmf * b * cosb * sinlam;
 140                         break;
 141                     case n_pole:
 142                     case s_pole:
 143                         if (q >= 0.) {
 144                             b = sqrt(q);
 145                             xy_x = b * sinlam;
 146                             xy_y = coslam * (this->m_proj_parm.mode == s_pole ? b : -b);
 147                         } else
 148                             xy_x = xy_y = 0.;
 149                         break;
 150                     }
 151                 }
 152
 153                 // INVERSE(e_inverse)  ellipsoid
 154                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 155                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 156                 {
 157                     T cCe, sCe, q, rho, ab=0.0;
 158
 159                     switch (this->m_proj_parm.mode) {
 160                     case equit:
 161                     case obliq:
 162                         xy_x /= this->m_proj_parm.dd;
 163                         xy_y *=  this->m_proj_parm.dd;
 164                         rho = boost::math::hypot(xy_x, xy_y);
 165                         if (rho < epsilon10) {
 166                             lp_lon = 0.;
 167                             lp_lat = par.phi0;
 168                             return;
 169                         }
 170                         sCe = 2. * asin(.5 * rho / this->m_proj_parm.rq);
 171                         cCe = cos(sCe);
 172                         sCe = sin(sCe);
 173                         xy_x *= sCe;
 174                         if (this->m_proj_parm.mode == obliq) {
 175                             ab = cCe * this->m_proj_parm.sinb1 + xy_y * sCe * this->m_proj_parm.cosb1 / rho;
 176                             xy_y = rho * this->m_proj_parm.cosb1 * cCe - xy_y * this->m_proj_parm.sinb1 * sCe;
 177                         } else {
 178                             ab = xy_y * sCe / rho;
 179                             xy_y = rho * cCe;
 180                         }
 181                         break;
 182                     case n_pole:
 183                         xy_y = -xy_y;
 184                         BOOST_FALLTHROUGH;
 185                     case s_pole:
 186                         q = (xy_x * xy_x + xy_y * xy_y);
 187                         if (q == 0.0) {
 188                             lp_lon = 0.;
 189                             lp_lat = par.phi0;
 190                             return;
 191                         }
 192                         ab = 1. - q / this->m_proj_parm.qp;
 193                         if (this->m_proj_parm.mode == s_pole)
 194                             ab = - ab;
 195                         break;
 196                     }
 197                     lp_lon = atan2(xy_x, xy_y);
 198                     lp_lat = pj_authlat(asin(ab), this->m_proj_parm.apa);
 199                 }
 200
 201                 static inline std::string get_name()
 202                 {
 203                     return "laea_ellipsoid";
 204                 }
 205
 206             };
 207
 208             template <typename T, typename Parameters>
 209             struct base_laea_spheroid
 210             {
 211                 par_laea<T> m_proj_parm;
 212
 213                 // FORWARD(s_forward)  spheroid
 214                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 215                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 216                 {
 217                     static const T fourth_pi = detail::fourth_pi<T>();
 218
 219                     T  coslam, cosphi, sinphi;
 220
 221                     sinphi = sin(lp_lat);
 222                     cosphi = cos(lp_lat);
 223                     coslam = cos(lp_lon);
 224                     switch (this->m_proj_parm.mode) {
 225                     case equit:
 226                         xy_y = 1. + cosphi * coslam;
 227                         goto oblcon;
 228                     case obliq:
 229                         xy_y = 1. + this->m_proj_parm.sinb1 * sinphi + this->m_proj_parm.cosb1 * cosphi * coslam;
 230                 oblcon:
 231                         if (xy_y <= epsilon10) {
 232                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 233                         }
 234                         xy_y = sqrt(2. / xy_y);
 235                         xy_x = xy_y * cosphi * sin(lp_lon);
 236                         xy_y *= this->m_proj_parm.mode == equit ? sinphi :
 237                            this->m_proj_parm.cosb1 * sinphi - this->m_proj_parm.sinb1 * cosphi * coslam;
 238                         break;
 239                     case n_pole:
 240                         coslam = -coslam;
 241                         BOOST_FALLTHROUGH;
 242                     case s_pole:
 243                         if (fabs(lp_lat + par.phi0) < epsilon10) {
 244                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 245                         }
 246                         xy_y = fourth_pi - lp_lat * .5;
 247                         xy_y = 2. * (this->m_proj_parm.mode == s_pole ? cos(xy_y) : sin(xy_y));
 248                         xy_x = xy_y * sin(lp_lon);
 249                         xy_y *= coslam;
 250                         break;
 251                     }
 252                 }
 253
 254                 // INVERSE(s_inverse)  spheroid
 255                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 256                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 257                 {
 258                     static const T half_pi = detail::half_pi<T>();
 259
 260                     T  cosz=0.0, rh, sinz=0.0;
 261
 262                     rh = boost::math::hypot(xy_x, xy_y);
 263                     if ((lp_lat = rh * .5 ) > 1.) {
 264                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 265                     }
 266                     lp_lat = 2. * asin(lp_lat);
 267                     if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
 268                         sinz = sin(lp_lat);
 269                         cosz = cos(lp_lat);
 270                     }
 271                     switch (this->m_proj_parm.mode) {
 272                     case equit:
 273                         lp_lat = fabs(rh) <= epsilon10 ? 0. : asin(xy_y * sinz / rh);
 274                         xy_x *= sinz;
 275                         xy_y = cosz * rh;
 276                         break;
 277                     case obliq:
 278                         lp_lat = fabs(rh) <= epsilon10 ? par.phi0 :
 279                            asin(cosz * this->m_proj_parm.sinb1 + xy_y * sinz * this->m_proj_parm.cosb1 / rh);
 280                         xy_x *= sinz * this->m_proj_parm.cosb1;
 281                         xy_y = (cosz - sin(lp_lat) * this->m_proj_parm.sinb1) * rh;
 282                         break;
 283                     case n_pole:
 284                         xy_y = -xy_y;
 285                         lp_lat = half_pi - lp_lat;
 286                         break;
 287                     case s_pole:
 288                         lp_lat -= half_pi;
 289                         break;
 290                     }
 291                     lp_lon = (xy_y == 0. && (this->m_proj_parm.mode == equit || this->m_proj_parm.mode == obliq)) ?
 292                         0. : atan2(xy_x, xy_y);
 293                 }
 294
 295                 static inline std::string get_name()
 296                 {
 297                     return "laea_spheroid";
 298                 }
 299
 300             };
 301
 302             // Lambert Azimuthal Equal Area
 303             template <typename Parameters, typename T>
 304             inline void setup_laea(Parameters& par, par_laea<T>& proj_parm)
 305             {
 306                 static const T half_pi = detail::half_pi<T>();
 307
 308                 T t;
 309
 310                 t = fabs(par.phi0);
 311                 if (fabs(t - half_pi) < epsilon10)
 312                     proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
 313                 else if (fabs(t) < epsilon10)
 314                     proj_parm.mode = equit;
 315                 else
 316                     proj_parm.mode = obliq;
 317                 if (par.es != 0.0) {
 318                     double sinphi;
 319
 320                     par.e = sqrt(par.es); // TODO : Isn't it already set?
 321                     proj_parm.qp = pj_qsfn(1., par.e, par.one_es);
 322                     proj_parm.mmf = .5 / (1. - par.es);
 323                     proj_parm.apa = pj_authset<T>(par.es);
 324                     switch (proj_parm.mode) {
 325                     case n_pole:
 326                     case s_pole:
 327                         proj_parm.dd = 1.;
 328                         break;
 329                     case equit:
 330                         proj_parm.dd = 1. / (proj_parm.rq = sqrt(.5 * proj_parm.qp));
 331                         proj_parm.xmf = 1.;
 332                         proj_parm.ymf = .5 * proj_parm.qp;
 333                         break;
 334                     case obliq:
 335                         proj_parm.rq = sqrt(.5 * proj_parm.qp);
 336                         sinphi = sin(par.phi0);
 337                         proj_parm.sinb1 = pj_qsfn(sinphi, par.e, par.one_es) / proj_parm.qp;
 338                         proj_parm.cosb1 = sqrt(1. - proj_parm.sinb1 * proj_parm.sinb1);
 339                         proj_parm.dd = cos(par.phi0) / (sqrt(1. - par.es * sinphi * sinphi) *
 340                            proj_parm.rq * proj_parm.cosb1);
 341                         proj_parm.ymf = (proj_parm.xmf = proj_parm.rq) / proj_parm.dd;
 342                         proj_parm.xmf *= proj_parm.dd;
 343                         break;
 344                     }
 345                 } else {
 346                     if (proj_parm.mode == obliq) {
 347                         proj_parm.sinb1 = sin(par.phi0);
 348                         proj_parm.cosb1 = cos(par.phi0);
 349                     }
 350                 }
 351             }
 352
 353     }} // namespace laea
 354     #endif // doxygen
 355
 356     /*!
 357         \brief Lambert Azimuthal Equal Area projection
 358         \ingroup projections
 359         \tparam Geographic latlong point type
 360         \tparam Cartesian xy point type
 361         \tparam Parameters parameter type
 362         \par Projection characteristics
 363          - Azimuthal
 364          - Spheroid
 365          - Ellipsoid
 366         \par Example
 367         \image html ex_laea.gif
 368     */
 369     template <typename T, typename Parameters>
 370     struct laea_ellipsoid : public detail::laea::base_laea_ellipsoid<T, Parameters>
 371     {
 372         template <typename Params>
 373         inline laea_ellipsoid(Params const& , Parameters & par)
 374         {
 375             detail::laea::setup_laea(par, this->m_proj_parm);
 376         }
 377     };
 378
 379     /*!
 380         \brief Lambert Azimuthal Equal Area projection
 381         \ingroup projections
 382         \tparam Geographic latlong point type
 383         \tparam Cartesian xy point type
 384         \tparam Parameters parameter type
 385         \par Projection characteristics
 386          - Azimuthal
 387          - Spheroid
 388          - Ellipsoid
 389         \par Example
 390         \image html ex_laea.gif
 391     */
 392     template <typename T, typename Parameters>
 393     struct laea_spheroid : public detail::laea::base_laea_spheroid<T, Parameters>
 394     {
 395         template <typename Params>
 396         inline laea_spheroid(Params const& , Parameters & par)
 397         {
 398             detail::laea::setup_laea(par, this->m_proj_parm);
 399         }
 400     };
 401
 402     #ifndef DOXYGEN_NO_DETAIL
 403     namespace detail
 404     {
 405
 406         // Static projection
 407         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_laea, laea_spheroid, laea_ellipsoid)
 408
 409         // Factory entry(s)
 410         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(laea_entry, laea_spheroid, laea_ellipsoid)
 411
 412         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(laea_init)
 413         {
 414             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(laea, laea_entry)
 415         }
 416
 417     } // namespace detail
 418     #endif // doxygen
 419
 420 } // namespace projections
 421
 422 }} // namespace boost::geometry
 423
 424 #endif // BOOST_GEOMETRY_PROJECTIONS_LAEA_HPP
 425