projections/proj/mod_ster.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_MOD_STER_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
0042
0043 #include <boost/geometry/util/math.hpp>
0044 #include <boost/math/special_functions/hypot.hpp>
0045
0046 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0047 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0048 #include <boost/geometry/srs/projections/impl/projects.hpp>
0049 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0050 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
0051 #include <boost/geometry/srs/projections/impl/pj_zpoly1.hpp>
0052
0053 namespace boost { namespace geometry
0054 {
0055
0056 namespace projections
0057 {
0058     #ifndef DOXYGEN_NO_DETAIL
0059     namespace detail { namespace mod_ster
0060     {
0061
0062             static const double epsilon = 1e-12;
0063
0064             template <typename T>
0065             struct par_mod_ster
0066             {
0067                 T              cchio, schio;
0068                 pj_complex<T>* zcoeff;
0069                 int            n;
0070             };
0071
0072             /* based upon Snyder and Linck, USGS-NMD */
0073
0074             template <typename T, typename Parameters>
0075             struct base_mod_ster_ellipsoid
0076             {
0077                 par_mod_ster<T> m_proj_parm;
0078
0079                 // FORWARD(e_forward)  ellipsoid
0080                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0081                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0082                 {
0083                     static const T half_pi = detail::half_pi<T>();
0084
0085                     T sinlon, coslon, esphi, chi, schi, cchi, s;
0086                     pj_complex<T> p;
0087
0088                     sinlon = sin(lp_lon);
0089                     coslon = cos(lp_lon);
0090                     esphi = par.e * sin(lp_lat);
0091                     chi = 2. * atan(tan((half_pi + lp_lat) * .5) *
0092                         math::pow((T(1) - esphi) / (T(1) + esphi), par.e * T(0.5))) - half_pi;
0093                     schi = sin(chi);
0094                     cchi = cos(chi);
0095                     s = 2. / (1. + this->m_proj_parm.schio * schi + this->m_proj_parm.cchio * cchi * coslon);
0096                     p.r = s * cchi * sinlon;
0097                     p.i = s * (this->m_proj_parm.cchio * schi - this->m_proj_parm.schio * cchi * coslon);
0098                     p = pj_zpoly1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n);
0099                     xy_x = p.r;
0100                     xy_y = p.i;
0101                 }
0102
0103                 // INVERSE(e_inverse)  ellipsoid
0104                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0105                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0106                 {
0107                     static const T half_pi = detail::half_pi<T>();
0108
0109                     int nn;
0110                     pj_complex<T> p, fxy, fpxy, dp;
0111                     T den, rh = 0, z, sinz = 0, cosz = 0, chi, phi = 0, dphi, esphi;
0112
0113                     p.r = xy_x;
0114                     p.i = xy_y;
0115                     for (nn = 20; nn ;--nn) {
0116                         fxy = pj_zpolyd1(p, this->m_proj_parm.zcoeff, this->m_proj_parm.n, &fpxy);
0117                         fxy.r -= xy_x;
0118                         fxy.i -= xy_y;
0119                         den = fpxy.r * fpxy.r + fpxy.i * fpxy.i;
0120                         dp.r = -(fxy.r * fpxy.r + fxy.i * fpxy.i) / den;
0121                         dp.i = -(fxy.i * fpxy.r - fxy.r * fpxy.i) / den;
0122                         p.r += dp.r;
0123                         p.i += dp.i;
0124                         if ((fabs(dp.r) + fabs(dp.i)) <= epsilon)
0125                             break;
0126                     }
0127                     if (nn) {
0128                         rh = boost::math::hypot(p.r, p.i);
0129                         z = 2. * atan(.5 * rh);
0130                         sinz = sin(z);
0131                         cosz = cos(z);
0132                         lp_lon = par.lam0;
0133                         if (fabs(rh) <= epsilon) {
0134                             /* if we end up here input coordinates were (0,0).
0135                              * pj_inv() adds P->lam0 to lp.lam, this way we are
0136                              * sure to get the correct offset */
0137                             lp_lon = 0.0;
0138                             lp_lat = par.phi0;
0139                             return;
0140                         }
0141                         chi = aasin(cosz * this->m_proj_parm.schio + p.i * sinz * this->m_proj_parm.cchio / rh);
0142                         phi = chi;
0143                         for (nn = 20; nn ;--nn) {
0144                             esphi = par.e * sin(phi);
0145                             dphi = 2. * atan(tan((half_pi + chi) * .5) *
0146                                 math::pow((T(1) + esphi) / (T(1) - esphi), par.e * T(0.5))) - half_pi - phi;
0147                             phi += dphi;
0148                             if (fabs(dphi) <= epsilon)
0149                                 break;
0150                         }
0151                     }
0152                     if (nn) {
0153                         lp_lat = phi;
0154                         lp_lon = atan2(p.r * sinz, rh * this->m_proj_parm.cchio * cosz - p.i *
0155                             this->m_proj_parm.schio * sinz);
0156                     } else
0157                         lp_lon = lp_lat = HUGE_VAL;
0158                 }
0159
0160                 static inline std::string get_name()
0161                 {
0162                     return "mod_ster_ellipsoid";
0163                 }
0164
0165             };
0166
0167             template <typename Parameters, typename T>
0168             inline void setup(Parameters& par, par_mod_ster<T>& proj_parm)  /* general initialization */
0169             {
0170                 static T const half_pi = detail::half_pi<T>();
0171
0172                 T esphi, chio;
0173
0174                 if (par.es != 0.0) {
0175                     esphi = par.e * sin(par.phi0);
0176                     chio = 2. * atan(tan((half_pi + par.phi0) * .5) *
0177                         math::pow((T(1) - esphi) / (T(1) + esphi), par.e * T(0.5))) - half_pi;
0178                 } else
0179                     chio = par.phi0;
0180                 proj_parm.schio = sin(chio);
0181                 proj_parm.cchio = cos(chio);
0182             }
0183
0184
0185             /* Miller Oblated Stereographic */
0186             template <typename Parameters, typename T>
0187             inline void setup_mil_os(Parameters& par, par_mod_ster<T>& proj_parm)
0188             {
0189                 static const T d2r = geometry::math::d2r<T>();
0190
0191                 static pj_complex<T> AB[] = {
0192                     {0.924500, 0.},
0193                     {0.,       0.},
0194                     {0.019430, 0.}
0195                 };
0196
0197                 proj_parm.n = 2;
0198                 par.lam0 = d2r * 20.;
0199                 par.phi0 = d2r * 18.;
0200                 proj_parm.zcoeff = AB;
0201                 par.es = 0.;
0202
0203                 setup(par, proj_parm);
0204             }
0205
0206             /* Lee Oblated Stereographic */
0207             template <typename Parameters, typename T>
0208             inline void setup_lee_os(Parameters& par, par_mod_ster<T>& proj_parm)
0209             {
0210                 static const T d2r = geometry::math::d2r<T>();
0211
0212                 static pj_complex<T> AB[] = {
0213                     { 0.721316,   0.},
0214                     { 0.,         0.},
0215                     {-0.0088162, -0.00617325}
0216                 };
0217
0218                 proj_parm.n = 2;
0219                 par.lam0 = d2r * -165.;
0220                 par.phi0 = d2r * -10.;
0221                 proj_parm.zcoeff = AB;
0222                 par.es = 0.;
0223
0224                 setup(par, proj_parm);
0225             }
0226
0227             // Mod. Stererographics of 48 U.S.
0228             template <typename Parameters, typename T>
0229             inline void setup_gs48(Parameters& par, par_mod_ster<T>& proj_parm)
0230             {
0231                 static const T d2r = geometry::math::d2r<T>();
0232
0233                 static pj_complex<T> AB[] = { /* 48 United States */
0234                     { 0.98879,  0.},
0235                     { 0.,       0.},
0236                     {-0.050909, 0.},
0237                     { 0.,       0.},
0238                     { 0.075528, 0.}
0239                 };
0240
0241                 proj_parm.n = 4;
0242                 par.lam0 = d2r * -96.;
0243                 par.phi0 = d2r * -39.;
0244                 proj_parm.zcoeff = AB;
0245                 par.es = 0.;
0246                 par.a = 6370997.;
0247
0248                 setup(par, proj_parm);
0249             }
0250
0251             // Mod. Stererographics of Alaska
0252             template <typename Parameters, typename T>
0253             inline void setup_alsk(Parameters& par, par_mod_ster<T>& proj_parm)
0254             {
0255                 static const T d2r = geometry::math::d2r<T>();
0256
0257                 static pj_complex<T> ABe[] = { /* Alaska ellipsoid */
0258                     { .9945303, 0.},
0259                     { .0052083, -.0027404},
0260                     { .0072721,  .0048181},
0261                     {-.0151089, -.1932526},
0262                     { .0642675, -.1381226},
0263                     { .3582802, -.2884586}
0264                 };
0265
0266                 static pj_complex<T> ABs[] = { /* Alaska sphere */
0267                     { .9972523, 0.},
0268                     { .0052513, -.0041175},
0269                     { .0074606,  .0048125},
0270                     {-.0153783, -.1968253},
0271                     { .0636871, -.1408027},
0272                     { .3660976, -.2937382}
0273                 };
0274
0275                 proj_parm.n = 5;
0276                 par.lam0 = d2r * -152.;
0277                 par.phi0 = d2r * 64.;
0278                 if (par.es != 0.0) { /* fixed ellipsoid/sphere */
0279                     proj_parm.zcoeff = ABe;
0280                     par.a = 6378206.4;
0281                     par.e = sqrt(par.es = 0.00676866);
0282                 } else {
0283                     proj_parm.zcoeff = ABs;
0284                     par.a = 6370997.;
0285                 }
0286
0287                 setup(par, proj_parm);
0288             }
0289
0290             // Mod. Stererographics of 50 U.S.
0291             template <typename Parameters, typename T>
0292             inline void setup_gs50(Parameters& par, par_mod_ster<T>& proj_parm)
0293             {
0294                 static const T d2r = geometry::math::d2r<T>();
0295
0296                 static pj_complex<T> ABe[] = { /* GS50 ellipsoid */
0297                     { .9827497, 0.},
0298                     { .0210669,  .0053804},
0299                     {-.1031415, -.0571664},
0300                     {-.0323337, -.0322847},
0301                     { .0502303,  .1211983},
0302                     { .0251805,  .0895678},
0303                     {-.0012315, -.1416121},
0304                     { .0072202, -.1317091},
0305                     {-.0194029,  .0759677},
0306                     {-.0210072,  .0834037}
0307                 };
0308                 static pj_complex<T> ABs[] = { /* GS50 sphere */
0309                     { .9842990, 0.},
0310                     { .0211642,  .0037608},
0311                     {-.1036018, -.0575102},
0312                     {-.0329095, -.0320119},
0313                     { .0499471,  .1223335},
0314                     { .0260460,  .0899805},
0315                     { .0007388, -.1435792},
0316                     { .0075848, -.1334108},
0317                     {-.0216473,  .0776645},
0318                     {-.0225161,  .0853673}
0319                 };
0320
0321                 proj_parm.n = 9;
0322                 par.lam0 = d2r * -120.;
0323                 par.phi0 = d2r * 45.;
0324                 if (par.es != 0.0) { /* fixed ellipsoid/sphere */
0325                     proj_parm.zcoeff = ABe;
0326                     par.a = 6378206.4;
0327                     par.e = sqrt(par.es = 0.00676866);
0328                 } else {
0329                     proj_parm.zcoeff = ABs;
0330                     par.a = 6370997.;
0331                 }
0332
0333                 setup(par, proj_parm);
0334             }
0335
0336     }} // namespace detail::mod_ster
0337     #endif // doxygen
0338
0339     /*!
0340         \brief Miller Oblated Stereographic projection
0341         \ingroup projections
0342         \tparam Geographic latlong point type
0343         \tparam Cartesian xy point type
0344         \tparam Parameters parameter type
0345         \par Projection characteristics
0346          - Azimuthal (mod)
0347         \par Example
0348         \image html ex_mil_os.gif
0349     */
0350     template <typename T, typename Parameters>
0351     struct mil_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
0352     {
0353         template <typename Params>
0354         inline mil_os_ellipsoid(Params const& , Parameters & par)
0355         {
0356             detail::mod_ster::setup_mil_os(par, this->m_proj_parm);
0357         }
0358     };
0359
0360     /*!
0361         \brief Lee Oblated Stereographic projection
0362         \ingroup projections
0363         \tparam Geographic latlong point type
0364         \tparam Cartesian xy point type
0365         \tparam Parameters parameter type
0366         \par Projection characteristics
0367          - Azimuthal (mod)
0368         \par Example
0369         \image html ex_lee_os.gif
0370     */
0371     template <typename T, typename Parameters>
0372     struct lee_os_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
0373     {
0374         template <typename Params>
0375         inline lee_os_ellipsoid(Params const& , Parameters & par)
0376         {
0377             detail::mod_ster::setup_lee_os(par, this->m_proj_parm);
0378         }
0379     };
0380
0381     /*!
0382         \brief Mod. Stererographics of 48 U.S. projection
0383         \ingroup projections
0384         \tparam Geographic latlong point type
0385         \tparam Cartesian xy point type
0386         \tparam Parameters parameter type
0387         \par Projection characteristics
0388          - Azimuthal (mod)
0389         \par Example
0390         \image html ex_gs48.gif
0391     */
0392     template <typename T, typename Parameters>
0393     struct gs48_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
0394     {
0395         template <typename Params>
0396         inline gs48_ellipsoid(Params const& , Parameters & par)
0397         {
0398             detail::mod_ster::setup_gs48(par, this->m_proj_parm);
0399         }
0400     };
0401
0402     /*!
0403         \brief Mod. Stererographics of Alaska projection
0404         \ingroup projections
0405         \tparam Geographic latlong point type
0406         \tparam Cartesian xy point type
0407         \tparam Parameters parameter type
0408         \par Projection characteristics
0409          - Azimuthal (mod)
0410         \par Example
0411         \image html ex_alsk.gif
0412     */
0413     template <typename T, typename Parameters>
0414     struct alsk_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
0415     {
0416         template <typename Params>
0417         inline alsk_ellipsoid(Params const& , Parameters & par)
0418         {
0419             detail::mod_ster::setup_alsk(par, this->m_proj_parm);
0420         }
0421     };
0422
0423     /*!
0424         \brief Mod. Stererographics of 50 U.S. projection
0425         \ingroup projections
0426         \tparam Geographic latlong point type
0427         \tparam Cartesian xy point type
0428         \tparam Parameters parameter type
0429         \par Projection characteristics
0430          - Azimuthal (mod)
0431         \par Example
0432         \image html ex_gs50.gif
0433     */
0434     template <typename T, typename Parameters>
0435     struct gs50_ellipsoid : public detail::mod_ster::base_mod_ster_ellipsoid<T, Parameters>
0436     {
0437         template <typename Params>
0438         inline gs50_ellipsoid(Params const& , Parameters & par)
0439         {
0440             detail::mod_ster::setup_gs50(par, this->m_proj_parm);
0441         }
0442     };
0443
0444     #ifndef DOXYGEN_NO_DETAIL
0445     namespace detail
0446     {
0447
0448         // Static projection
0449         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_mil_os, mil_os_ellipsoid)
0450         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lee_os, lee_os_ellipsoid)
0451         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_gs48, gs48_ellipsoid)
0452         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_alsk, alsk_ellipsoid)
0453         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_gs50, gs50_ellipsoid)
0454
0455         // Factory entry(s)
0456         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(mil_os_entry, mil_os_ellipsoid)
0457         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lee_os_entry, lee_os_ellipsoid)
0458         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(gs48_entry, gs48_ellipsoid)
0459         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(alsk_entry, alsk_ellipsoid)
0460         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(gs50_entry, gs50_ellipsoid)
0461
0462         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(mod_ster_init)
0463         {
0464             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(mil_os, mil_os_entry)
0465             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lee_os, lee_os_entry)
0466             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(gs48, gs48_entry)
0467             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(alsk, alsk_entry)
0468             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(gs50, gs50_entry)
0469         }
0470
0471     } // namespace detail
0472     #endif // doxygen
0473
0474 } // namespace projections
0475
0476 }} // namespace boost::geometry
0477
0478 #endif // BOOST_GEOMETRY_PROJECTIONS_MOD_STER_HPP
0479