projections/proj/laea.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Permission is hereby granted, free of charge, to any person obtaining a
0023 // copy of this software and associated documentation files (the "Software"),
0024 // to deal in the Software without restriction, including without limitation
0025 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0026 // and/or sell copies of the Software, and to permit persons to whom the
0027 // Software is furnished to do so, subject to the following conditions:
0028
0029 // The above copyright notice and this permission notice shall be included
0030 // in all copies or substantial portions of the Software.
0031
0032 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0033 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0034 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0035 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0036 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0037 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0038 // DEALINGS IN THE SOFTWARE.
0039
0040 #ifndef BOOST_GEOMETRY_PROJECTIONS_LAEA_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_LAEA_HPP
0042
0043 #include <boost/config.hpp>
0044 #include <boost/geometry/util/math.hpp>
0045 #include <boost/math/special_functions/hypot.hpp>
0046
0047 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0048 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0049 #include <boost/geometry/srs/projections/impl/projects.hpp>
0050 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0051 #include <boost/geometry/srs/projections/impl/pj_auth.hpp>
0052 #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
0053
0054 namespace boost { namespace geometry
0055 {
0056
0057 namespace projections
0058 {
0059     #ifndef DOXYGEN_NO_DETAIL
0060     namespace detail { namespace laea
0061     {
0062             static const double epsilon10 = 1.e-10;
0063
0064             enum mode_type {
0065                 n_pole = 0,
0066                 s_pole = 1,
0067                 equit  = 2,
0068                 obliq  = 3
0069             };
0070
0071             template <typename T>
0072             struct par_laea
0073             {
0074                 T   sinb1;
0075                 T   cosb1;
0076                 T   xmf;
0077                 T   ymf;
0078                 T   mmf;
0079                 T   qp;
0080                 T   dd;
0081                 T   rq;
0082                 detail::apa<T> apa;
0083                 mode_type mode;
0084             };
0085
0086             template <typename T, typename Parameters>
0087             struct base_laea_ellipsoid
0088             {
0089                 par_laea<T> m_proj_parm;
0090
0091                 // FORWARD(e_forward)  ellipsoid
0092                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0093                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0094                 {
0095                     static const T half_pi = detail::half_pi<T>();
0096
0097                     T coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
0098
0099                     coslam = cos(lp_lon);
0100                     sinlam = sin(lp_lon);
0101                     sinphi = sin(lp_lat);
0102                     q = pj_qsfn(sinphi, par.e, par.one_es);
0103
0104                     if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
0105                         sinb = q / this->m_proj_parm.qp;
0106                         cosb = sqrt(1. - sinb * sinb);
0107                     }
0108
0109                     switch (this->m_proj_parm.mode) {
0110                     case obliq:
0111                         b = 1. + this->m_proj_parm.sinb1 * sinb + this->m_proj_parm.cosb1 * cosb * coslam;
0112                         break;
0113                     case equit:
0114                         b = 1. + cosb * coslam;
0115                         break;
0116                     case n_pole:
0117                         b = half_pi + lp_lat;
0118                         q = this->m_proj_parm.qp - q;
0119                         break;
0120                     case s_pole:
0121                         b = lp_lat - half_pi;
0122                         q = this->m_proj_parm.qp + q;
0123                         break;
0124                     }
0125                     if (fabs(b) < epsilon10) {
0126                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0127                     }
0128
0129                     switch (this->m_proj_parm.mode) {
0130                     case obliq:
0131                         b = sqrt(2. / b);
0132                         xy_y = this->m_proj_parm.ymf * b * (this->m_proj_parm.cosb1 * sinb - this->m_proj_parm.sinb1 * cosb * coslam);
0133                         goto eqcon;
0134                         break;
0135                     case equit:
0136                         b = sqrt(2. / (1. + cosb * coslam));
0137                         xy_y = b * sinb * this->m_proj_parm.ymf;
0138                 eqcon:
0139                         xy_x = this->m_proj_parm.xmf * b * cosb * sinlam;
0140                         break;
0141                     case n_pole:
0142                     case s_pole:
0143                         if (q >= 0.) {
0144                             b = sqrt(q);
0145                             xy_x = b * sinlam;
0146                             xy_y = coslam * (this->m_proj_parm.mode == s_pole ? b : -b);
0147                         } else
0148                             xy_x = xy_y = 0.;
0149                         break;
0150                     }
0151                 }
0152
0153                 // INVERSE(e_inverse)  ellipsoid
0154                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0155                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0156                 {
0157                     T cCe, sCe, q, rho, ab=0.0;
0158
0159                     switch (this->m_proj_parm.mode) {
0160                     case equit:
0161                     case obliq:
0162                         xy_x /= this->m_proj_parm.dd;
0163                         xy_y *=  this->m_proj_parm.dd;
0164                         rho = boost::math::hypot(xy_x, xy_y);
0165                         if (rho < epsilon10) {
0166                             lp_lon = 0.;
0167                             lp_lat = par.phi0;
0168                             return;
0169                         }
0170                         sCe = 2. * asin(.5 * rho / this->m_proj_parm.rq);
0171                         cCe = cos(sCe);
0172                         sCe = sin(sCe);
0173                         xy_x *= sCe;
0174                         if (this->m_proj_parm.mode == obliq) {
0175                             ab = cCe * this->m_proj_parm.sinb1 + xy_y * sCe * this->m_proj_parm.cosb1 / rho;
0176                             xy_y = rho * this->m_proj_parm.cosb1 * cCe - xy_y * this->m_proj_parm.sinb1 * sCe;
0177                         } else {
0178                             ab = xy_y * sCe / rho;
0179                             xy_y = rho * cCe;
0180                         }
0181                         break;
0182                     case n_pole:
0183                         xy_y = -xy_y;
0184                         BOOST_FALLTHROUGH;
0185                     case s_pole:
0186                         q = (xy_x * xy_x + xy_y * xy_y);
0187                         if (q == 0.0) {
0188                             lp_lon = 0.;
0189                             lp_lat = par.phi0;
0190                             return;
0191                         }
0192                         ab = 1. - q / this->m_proj_parm.qp;
0193                         if (this->m_proj_parm.mode == s_pole)
0194                             ab = - ab;
0195                         break;
0196                     }
0197                     lp_lon = atan2(xy_x, xy_y);
0198                     lp_lat = pj_authlat(asin(ab), this->m_proj_parm.apa);
0199                 }
0200
0201                 static inline std::string get_name()
0202                 {
0203                     return "laea_ellipsoid";
0204                 }
0205
0206             };
0207
0208             template <typename T, typename Parameters>
0209             struct base_laea_spheroid
0210             {
0211                 par_laea<T> m_proj_parm;
0212
0213                 // FORWARD(s_forward)  spheroid
0214                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0215                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0216                 {
0217                     static const T fourth_pi = detail::fourth_pi<T>();
0218
0219                     T  coslam, cosphi, sinphi;
0220
0221                     sinphi = sin(lp_lat);
0222                     cosphi = cos(lp_lat);
0223                     coslam = cos(lp_lon);
0224                     switch (this->m_proj_parm.mode) {
0225                     case equit:
0226                         xy_y = 1. + cosphi * coslam;
0227                         goto oblcon;
0228                     case obliq:
0229                         xy_y = 1. + this->m_proj_parm.sinb1 * sinphi + this->m_proj_parm.cosb1 * cosphi * coslam;
0230                 oblcon:
0231                         if (xy_y <= epsilon10) {
0232                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0233                         }
0234                         xy_y = sqrt(2. / xy_y);
0235                         xy_x = xy_y * cosphi * sin(lp_lon);
0236                         xy_y *= this->m_proj_parm.mode == equit ? sinphi :
0237                            this->m_proj_parm.cosb1 * sinphi - this->m_proj_parm.sinb1 * cosphi * coslam;
0238                         break;
0239                     case n_pole:
0240                         coslam = -coslam;
0241                         BOOST_FALLTHROUGH;
0242                     case s_pole:
0243                         if (fabs(lp_lat + par.phi0) < epsilon10) {
0244                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0245                         }
0246                         xy_y = fourth_pi - lp_lat * .5;
0247                         xy_y = 2. * (this->m_proj_parm.mode == s_pole ? cos(xy_y) : sin(xy_y));
0248                         xy_x = xy_y * sin(lp_lon);
0249                         xy_y *= coslam;
0250                         break;
0251                     }
0252                 }
0253
0254                 // INVERSE(s_inverse)  spheroid
0255                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0256                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0257                 {
0258                     static const T half_pi = detail::half_pi<T>();
0259
0260                     T  cosz=0.0, rh, sinz=0.0;
0261
0262                     rh = boost::math::hypot(xy_x, xy_y);
0263                     if ((lp_lat = rh * .5 ) > 1.) {
0264                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0265                     }
0266                     lp_lat = 2. * asin(lp_lat);
0267                     if (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit) {
0268                         sinz = sin(lp_lat);
0269                         cosz = cos(lp_lat);
0270                     }
0271                     switch (this->m_proj_parm.mode) {
0272                     case equit:
0273                         lp_lat = fabs(rh) <= epsilon10 ? 0. : asin(xy_y * sinz / rh);
0274                         xy_x *= sinz;
0275                         xy_y = cosz * rh;
0276                         break;
0277                     case obliq:
0278                         lp_lat = fabs(rh) <= epsilon10 ? par.phi0 :
0279                            asin(cosz * this->m_proj_parm.sinb1 + xy_y * sinz * this->m_proj_parm.cosb1 / rh);
0280                         xy_x *= sinz * this->m_proj_parm.cosb1;
0281                         xy_y = (cosz - sin(lp_lat) * this->m_proj_parm.sinb1) * rh;
0282                         break;
0283                     case n_pole:
0284                         xy_y = -xy_y;
0285                         lp_lat = half_pi - lp_lat;
0286                         break;
0287                     case s_pole:
0288                         lp_lat -= half_pi;
0289                         break;
0290                     }
0291                     lp_lon = (xy_y == 0. && (this->m_proj_parm.mode == equit || this->m_proj_parm.mode == obliq)) ?
0292                         0. : atan2(xy_x, xy_y);
0293                 }
0294
0295                 static inline std::string get_name()
0296                 {
0297                     return "laea_spheroid";
0298                 }
0299
0300             };
0301
0302             // Lambert Azimuthal Equal Area
0303             template <typename Parameters, typename T>
0304             inline void setup_laea(Parameters& par, par_laea<T>& proj_parm)
0305             {
0306                 static const T half_pi = detail::half_pi<T>();
0307
0308                 T t;
0309
0310                 t = fabs(par.phi0);
0311                 if (fabs(t - half_pi) < epsilon10)
0312                     proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
0313                 else if (fabs(t) < epsilon10)
0314                     proj_parm.mode = equit;
0315                 else
0316                     proj_parm.mode = obliq;
0317                 if (par.es != 0.0) {
0318                     double sinphi;
0319
0320                     par.e = sqrt(par.es); // TODO : Isn't it already set?
0321                     proj_parm.qp = pj_qsfn(1., par.e, par.one_es);
0322                     proj_parm.mmf = .5 / (1. - par.es);
0323                     proj_parm.apa = pj_authset<T>(par.es);
0324                     switch (proj_parm.mode) {
0325                     case n_pole:
0326                     case s_pole:
0327                         proj_parm.dd = 1.;
0328                         break;
0329                     case equit:
0330                         proj_parm.dd = 1. / (proj_parm.rq = sqrt(.5 * proj_parm.qp));
0331                         proj_parm.xmf = 1.;
0332                         proj_parm.ymf = .5 * proj_parm.qp;
0333                         break;
0334                     case obliq:
0335                         proj_parm.rq = sqrt(.5 * proj_parm.qp);
0336                         sinphi = sin(par.phi0);
0337                         proj_parm.sinb1 = pj_qsfn(sinphi, par.e, par.one_es) / proj_parm.qp;
0338                         proj_parm.cosb1 = sqrt(1. - proj_parm.sinb1 * proj_parm.sinb1);
0339                         proj_parm.dd = cos(par.phi0) / (sqrt(1. - par.es * sinphi * sinphi) *
0340                            proj_parm.rq * proj_parm.cosb1);
0341                         proj_parm.ymf = (proj_parm.xmf = proj_parm.rq) / proj_parm.dd;
0342                         proj_parm.xmf *= proj_parm.dd;
0343                         break;
0344                     }
0345                 } else {
0346                     if (proj_parm.mode == obliq) {
0347                         proj_parm.sinb1 = sin(par.phi0);
0348                         proj_parm.cosb1 = cos(par.phi0);
0349                     }
0350                 }
0351             }
0352
0353     }} // namespace laea
0354     #endif // doxygen
0355
0356     /*!
0357         \brief Lambert Azimuthal Equal Area projection
0358         \ingroup projections
0359         \tparam Geographic latlong point type
0360         \tparam Cartesian xy point type
0361         \tparam Parameters parameter type
0362         \par Projection characteristics
0363          - Azimuthal
0364          - Spheroid
0365          - Ellipsoid
0366         \par Example
0367         \image html ex_laea.gif
0368     */
0369     template <typename T, typename Parameters>
0370     struct laea_ellipsoid : public detail::laea::base_laea_ellipsoid<T, Parameters>
0371     {
0372         template <typename Params>
0373         inline laea_ellipsoid(Params const& , Parameters & par)
0374         {
0375             detail::laea::setup_laea(par, this->m_proj_parm);
0376         }
0377     };
0378
0379     /*!
0380         \brief Lambert Azimuthal Equal Area projection
0381         \ingroup projections
0382         \tparam Geographic latlong point type
0383         \tparam Cartesian xy point type
0384         \tparam Parameters parameter type
0385         \par Projection characteristics
0386          - Azimuthal
0387          - Spheroid
0388          - Ellipsoid
0389         \par Example
0390         \image html ex_laea.gif
0391     */
0392     template <typename T, typename Parameters>
0393     struct laea_spheroid : public detail::laea::base_laea_spheroid<T, Parameters>
0394     {
0395         template <typename Params>
0396         inline laea_spheroid(Params const& , Parameters & par)
0397         {
0398             detail::laea::setup_laea(par, this->m_proj_parm);
0399         }
0400     };
0401
0402     #ifndef DOXYGEN_NO_DETAIL
0403     namespace detail
0404     {
0405
0406         // Static projection
0407         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_laea, laea_spheroid, laea_ellipsoid)
0408
0409         // Factory entry(s)
0410         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(laea_entry, laea_spheroid, laea_ellipsoid)
0411
0412         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(laea_init)
0413         {
0414             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(laea, laea_entry)
0415         }
0416
0417     } // namespace detail
0418     #endif // doxygen
0419
0420 } // namespace projections
0421
0422 }} // namespace boost::geometry
0423
0424 #endif // BOOST_GEOMETRY_PROJECTIONS_LAEA_HPP
0425