projections/proj/lsat.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_LSAT_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_LSAT_HPP
0042
0043 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
0044 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0045 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0046 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0047 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0048 #include <boost/geometry/srs/projections/impl/projects.hpp>
0049
0050 #include <boost/geometry/util/math.hpp>
0051
0052 namespace boost { namespace geometry
0053 {
0054
0055 namespace projections
0056 {
0057     #ifndef DOXYGEN_NO_DETAIL
0058     namespace detail { namespace lsat
0059     {
0060             static const double tolerance = 1e-7;
0061
0062             template <typename T>
0063             struct par_lsat
0064             {
0065                 T a2, a4, b, c1, c3;
0066                 T q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
0067             };
0068
0069             /* based upon Snyder and Linck, USGS-NMD */
0070             template <typename T>
0071             inline void
0072             seraz0(T lam, T const& mult, par_lsat<T>& proj_parm)
0073             {
0074                 T sdsq, h, s, fc, sd, sq, d__1 = 0;
0075
0076                 lam *= geometry::math::d2r<T>();
0077                 sd = sin(lam);
0078                 sdsq = sd * sd;
0079                 s = proj_parm.p22 * proj_parm.sa * cos(lam) * sqrt((1. + proj_parm.t * sdsq)
0080                     / ((1. + proj_parm.w * sdsq) * (1. + proj_parm.q * sdsq)));
0081
0082                 d__1 = 1. + proj_parm.q * sdsq;
0083                 h = sqrt((1. + proj_parm.q * sdsq) / (1. + proj_parm.w * sdsq)) * ((1. + proj_parm.w * sdsq)
0084                     / (d__1 * d__1) - proj_parm.p22 * proj_parm.ca);
0085
0086                 sq = sqrt(proj_parm.xj * proj_parm.xj + s * s);
0087                 fc = mult * (h * proj_parm.xj - s * s) / sq;
0088                 proj_parm.b += fc;
0089                 proj_parm.a2 += fc * cos(lam + lam);
0090                 proj_parm.a4 += fc * cos(lam * 4.);
0091                 fc = mult * s * (h + proj_parm.xj) / sq;
0092                 proj_parm.c1 += fc * cos(lam);
0093                 proj_parm.c3 += fc * cos(lam * 3.);
0094             }
0095
0096             template <typename T, typename Parameters>
0097             struct base_lsat_ellipsoid
0098             {
0099                 par_lsat<T> m_proj_parm;
0100
0101                 // FORWARD(e_forward)  ellipsoid
0102                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0103                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
0104                 {
0105                     static const T fourth_pi = detail::fourth_pi<T>();
0106                     static const T half_pi = detail::half_pi<T>();
0107                     static const T one_and_half_pi = detail::one_and_half_pi<T>();
0108                     static const T two_and_half_pi = detail::two_and_half_pi<T>();
0109
0110                     int l, nn;
0111                     T lamt = 0.0, xlam, sdsq, c, d, s, lamdp = 0.0, phidp, lampp, tanph;
0112                     T lamtp, cl, sd, sp, sav, tanphi;
0113
0114                     if (lp_lat > half_pi)
0115                         lp_lat = half_pi;
0116                     else if (lp_lat < -half_pi)
0117                         lp_lat = -half_pi;
0118
0119                     if (lp_lat >= 0. )
0120                         lampp = half_pi;
0121                     else
0122                         lampp = one_and_half_pi;
0123                     tanphi = tan(lp_lat);
0124                     for (nn = 0;;) {
0125                         T fac;
0126                         sav = lampp;
0127                         lamtp = lp_lon + this->m_proj_parm.p22 * lampp;
0128                         cl = cos(lamtp);
0129                         if (fabs(cl) < tolerance)
0130                             lamtp -= tolerance;
0131                         if( cl < 0 )
0132                             fac = lampp + sin(lampp) * half_pi;
0133                         else
0134                             fac = lampp - sin(lampp) * half_pi;
0135                         for (l = 50; l; --l) {
0136                             lamt = lp_lon + this->m_proj_parm.p22 * sav;
0137                             c = cos(lamt);
0138                             if (fabs(c) < tolerance)
0139                                 lamt -= tolerance;
0140                             xlam = (par.one_es * tanphi * this->m_proj_parm.sa + sin(lamt) * this->m_proj_parm.ca) / c;
0141                             lamdp = atan(xlam) + fac;
0142                             if (fabs(fabs(sav) - fabs(lamdp)) < tolerance)
0143                                 break;
0144                             sav = lamdp;
0145                         }
0146                         if (!l || ++nn >= 3 || (lamdp > this->m_proj_parm.rlm && lamdp < this->m_proj_parm.rlm2))
0147                             break;
0148                         if (lamdp <= this->m_proj_parm.rlm)
0149                             lampp = two_and_half_pi;
0150                         else if (lamdp >= this->m_proj_parm.rlm2)
0151                             lampp = half_pi;
0152                     }
0153                     if (l) {
0154                         sp = sin(lp_lat);
0155                         phidp = aasin((par.one_es * this->m_proj_parm.ca * sp - this->m_proj_parm.sa * cos(lp_lat) *
0156                             sin(lamt)) / sqrt(1. - par.es * sp * sp));
0157                         tanph = log(tan(fourth_pi + .5 * phidp));
0158                         sd = sin(lamdp);
0159                         sdsq = sd * sd;
0160                         s = this->m_proj_parm.p22 * this->m_proj_parm.sa * cos(lamdp) * sqrt((1. + this->m_proj_parm.t * sdsq)
0161                              / ((1. + this->m_proj_parm.w * sdsq) * (1. + this->m_proj_parm.q * sdsq)));
0162                         d = sqrt(this->m_proj_parm.xj * this->m_proj_parm.xj + s * s);
0163                         xy_x = this->m_proj_parm.b * lamdp + this->m_proj_parm.a2 * sin(2. * lamdp) + this->m_proj_parm.a4 *
0164                             sin(lamdp * 4.) - tanph * s / d;
0165                         xy_y = this->m_proj_parm.c1 * sd + this->m_proj_parm.c3 * sin(lamdp * 3.) + tanph * this->m_proj_parm.xj / d;
0166                     } else
0167                         xy_x = xy_y = HUGE_VAL;
0168                 }
0169
0170                 // INVERSE(e_inverse)  ellipsoid
0171                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0172                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0173                 {
0174                     static const T fourth_pi = detail::fourth_pi<T>();
0175                     static const T half_pi = detail::half_pi<T>();
0176
0177                     int nn;
0178                     T lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
0179
0180                     lamdp = xy_x / this->m_proj_parm.b;
0181                     nn = 50;
0182                     do {
0183                         sav = lamdp;
0184                         sd = sin(lamdp);
0185                         sdsq = sd * sd;
0186                         s = this->m_proj_parm.p22 * this->m_proj_parm.sa * cos(lamdp) * sqrt((1. + this->m_proj_parm.t * sdsq)
0187                              / ((1. + this->m_proj_parm.w * sdsq) * (1. + this->m_proj_parm.q * sdsq)));
0188                         lamdp = xy_x + xy_y * s / this->m_proj_parm.xj - this->m_proj_parm.a2 * sin(
0189                             2. * lamdp) - this->m_proj_parm.a4 * sin(lamdp * 4.) - s / this->m_proj_parm.xj * (
0190                             this->m_proj_parm.c1 * sin(lamdp) + this->m_proj_parm.c3 * sin(lamdp * 3.));
0191                         lamdp /= this->m_proj_parm.b;
0192                     } while (fabs(lamdp - sav) >= tolerance && --nn);
0193                     sl = sin(lamdp);
0194                     fac = exp(sqrt(1. + s * s / this->m_proj_parm.xj / this->m_proj_parm.xj) * (xy_y -
0195                         this->m_proj_parm.c1 * sl - this->m_proj_parm.c3 * sin(lamdp * 3.)));
0196                     phidp = 2. * (atan(fac) - fourth_pi);
0197                     dd = sl * sl;
0198                     if (fabs(cos(lamdp)) < tolerance)
0199                         lamdp -= tolerance;
0200                     spp = sin(phidp);
0201                     sppsq = spp * spp;
0202                     lamt = atan(((1. - sppsq * par.rone_es) * tan(lamdp) *
0203                         this->m_proj_parm.ca - spp * this->m_proj_parm.sa * sqrt((1. + this->m_proj_parm.q * dd) * (
0204                         1. - sppsq) - sppsq * this->m_proj_parm.u) / cos(lamdp)) / (1. - sppsq
0205                         * (1. + this->m_proj_parm.u)));
0206                     sl = lamt >= 0. ? 1. : -1.;
0207                     scl = cos(lamdp) >= 0. ? 1. : -1;
0208                     lamt -= half_pi * (1. - scl) * sl;
0209                     lp_lon = lamt - this->m_proj_parm.p22 * lamdp;
0210                     if (fabs(this->m_proj_parm.sa) < tolerance)
0211                         lp_lat = aasin(spp / sqrt(par.one_es * par.one_es + par.es * sppsq));
0212                     else
0213                         lp_lat = atan((tan(lamdp) * cos(lamt) - this->m_proj_parm.ca * sin(lamt)) /
0214                             (par.one_es * this->m_proj_parm.sa));
0215                 }
0216
0217                 static inline std::string get_name()
0218                 {
0219                     return "lsat_ellipsoid";
0220                 }
0221
0222             };
0223
0224             // Space oblique for LANDSAT
0225             template <typename Params, typename Parameters, typename T>
0226             inline void setup_lsat(Params const& params, Parameters& par, par_lsat<T>& proj_parm)
0227             {
0228                 static T const d2r = geometry::math::d2r<T>();
0229                 static T const pi = detail::pi<T>();
0230                 static T const two_pi = detail::two_pi<T>();
0231
0232                 int land, path;
0233                 T lam, alf, esc, ess;
0234
0235                 land = pj_get_param_i<srs::spar::lsat>(params, "lsat", srs::dpar::lsat);
0236                 if (land <= 0 || land > 5)
0237                     BOOST_THROW_EXCEPTION( projection_exception(error_lsat_not_in_range) );
0238
0239                 path = pj_get_param_i<srs::spar::path>(params, "path", srs::dpar::path);
0240                 if (path <= 0 || path > (land <= 3 ? 251 : 233))
0241                     BOOST_THROW_EXCEPTION( projection_exception(error_path_not_in_range) );
0242
0243                 if (land <= 3) {
0244                     par.lam0 = d2r * 128.87 - two_pi / 251. * path;
0245                     proj_parm.p22 = 103.2669323;
0246                     alf = d2r * 99.092;
0247                 } else {
0248                     par.lam0 = d2r * 129.3 - two_pi / 233. * path;
0249                     proj_parm.p22 = 98.8841202;
0250                     alf = d2r * 98.2;
0251                 }
0252                 proj_parm.p22 /= 1440.;
0253                 proj_parm.sa = sin(alf);
0254                 proj_parm.ca = cos(alf);
0255                 if (fabs(proj_parm.ca) < 1e-9)
0256                     proj_parm.ca = 1e-9;
0257                 esc = par.es * proj_parm.ca * proj_parm.ca;
0258                 ess = par.es * proj_parm.sa * proj_parm.sa;
0259                 proj_parm.w = (1. - esc) * par.rone_es;
0260                 proj_parm.w = proj_parm.w * proj_parm.w - 1.;
0261                 proj_parm.q = ess * par.rone_es;
0262                 proj_parm.t = ess * (2. - par.es) * par.rone_es * par.rone_es;
0263                 proj_parm.u = esc * par.rone_es;
0264                 proj_parm.xj = par.one_es * par.one_es * par.one_es;
0265                 proj_parm.rlm = pi * (1. / 248. + .5161290322580645);
0266                 proj_parm.rlm2 = proj_parm.rlm + two_pi;
0267                 proj_parm.a2 = proj_parm.a4 = proj_parm.b = proj_parm.c1 = proj_parm.c3 = 0.;
0268                 seraz0(0., 1., proj_parm);
0269                 for (lam = 9.; lam <= 81.0001; lam += 18.)
0270                     seraz0(lam, 4., proj_parm);
0271                 for (lam = 18; lam <= 72.0001; lam += 18.)
0272                     seraz0(lam, 2., proj_parm);
0273                 seraz0(90., 1., proj_parm);
0274                 proj_parm.a2 /= 30.;
0275                 proj_parm.a4 /= 60.;
0276                 proj_parm.b /= 30.;
0277                 proj_parm.c1 /= 15.;
0278                 proj_parm.c3 /= 45.;
0279             }
0280
0281     }} // namespace detail::lsat
0282     #endif // doxygen
0283
0284     /*!
0285         \brief Space oblique for LANDSAT projection
0286         \ingroup projections
0287         \tparam Geographic latlong point type
0288         \tparam Cartesian xy point type
0289         \tparam Parameters parameter type
0290         \par Projection characteristics
0291          - Cylindrical
0292          - Spheroid
0293          - Ellipsoid
0294         \par Projection parameters
0295          - lsat (integer)
0296          - path (integer)
0297         \par Example
0298         \image html ex_lsat.gif
0299     */
0300     template <typename T, typename Parameters>
0301     struct lsat_ellipsoid : public detail::lsat::base_lsat_ellipsoid<T, Parameters>
0302     {
0303         template <typename Params>
0304         inline lsat_ellipsoid(Params const& params, Parameters & par)
0305         {
0306             detail::lsat::setup_lsat(params, par, this->m_proj_parm);
0307         }
0308     };
0309
0310     #ifndef DOXYGEN_NO_DETAIL
0311     namespace detail
0312     {
0313
0314         // Static projection
0315         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lsat, lsat_ellipsoid)
0316
0317         // Factory entry(s)
0318         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lsat_entry, lsat_ellipsoid)
0319
0320         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(lsat_init)
0321         {
0322             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lsat, lsat_entry)
0323         }
0324
0325     } // namespace detail
0326     #endif // doxygen
0327
0328 } // namespace projections
0329
0330 }} // namespace boost::geometry
0331
0332 #endif // BOOST_GEOMETRY_PROJECTIONS_LSAT_HPP
0333