boost/geometry/srs/projections/proj/lsat.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_LSAT_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_LSAT_HPP
  42
  43 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
  44 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  45 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  46 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  47 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  48 #include <boost/geometry/srs/projections/impl/projects.hpp>
  49
  50 #include <boost/geometry/util/math.hpp>
  51
  52 namespace boost { namespace geometry
  53 {
  54
  55 namespace projections
  56 {
  57     #ifndef DOXYGEN_NO_DETAIL
  58     namespace detail { namespace lsat
  59     {
  60             static const double tolerance = 1e-7;
  61
  62             template <typename T>
  63             struct par_lsat
  64             {
  65                 T a2, a4, b, c1, c3;
  66                 T q, t, u, w, p22, sa, ca, xj, rlm, rlm2;
  67             };
  68
  69             /* based upon Snyder and Linck, USGS-NMD */
  70             template <typename T>
  71             inline void
  72             seraz0(T lam, T const& mult, par_lsat<T>& proj_parm)
  73             {
  74                 T sdsq, h, s, fc, sd, sq, d__1 = 0;
  75
  76                 lam *= geometry::math::d2r<T>();
  77                 sd = sin(lam);
  78                 sdsq = sd * sd;
  79                 s = proj_parm.p22 * proj_parm.sa * cos(lam) * sqrt((1. + proj_parm.t * sdsq)
  80                     / ((1. + proj_parm.w * sdsq) * (1. + proj_parm.q * sdsq)));
  81
  82                 d__1 = 1. + proj_parm.q * sdsq;
  83                 h = sqrt((1. + proj_parm.q * sdsq) / (1. + proj_parm.w * sdsq)) * ((1. + proj_parm.w * sdsq)
  84                     / (d__1 * d__1) - proj_parm.p22 * proj_parm.ca);
  85
  86                 sq = sqrt(proj_parm.xj * proj_parm.xj + s * s);
  87                 fc = mult * (h * proj_parm.xj - s * s) / sq;
  88                 proj_parm.b += fc;
  89                 proj_parm.a2 += fc * cos(lam + lam);
  90                 proj_parm.a4 += fc * cos(lam * 4.);
  91                 fc = mult * s * (h + proj_parm.xj) / sq;
  92                 proj_parm.c1 += fc * cos(lam);
  93                 proj_parm.c3 += fc * cos(lam * 3.);
  94             }
  95
  96             template <typename T, typename Parameters>
  97             struct base_lsat_ellipsoid
  98             {
  99                 par_lsat<T> m_proj_parm;
 100
 101                 // FORWARD(e_forward)  ellipsoid
 102                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 103                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
 104                 {
 105                     static const T fourth_pi = detail::fourth_pi<T>();
 106                     static const T half_pi = detail::half_pi<T>();
 107                     static const T one_and_half_pi = detail::one_and_half_pi<T>();
 108                     static const T two_and_half_pi = detail::two_and_half_pi<T>();
 109
 110                     int l, nn;
 111                     T lamt = 0.0, xlam, sdsq, c, d, s, lamdp = 0.0, phidp, lampp, tanph;
 112                     T lamtp, cl, sd, sp, sav, tanphi;
 113
 114                     if (lp_lat > half_pi)
 115                         lp_lat = half_pi;
 116                     else if (lp_lat < -half_pi)
 117                         lp_lat = -half_pi;
 118
 119                     if (lp_lat >= 0. )
 120                         lampp = half_pi;
 121                     else
 122                         lampp = one_and_half_pi;
 123                     tanphi = tan(lp_lat);
 124                     for (nn = 0;;) {
 125                         T fac;
 126                         sav = lampp;
 127                         lamtp = lp_lon + this->m_proj_parm.p22 * lampp;
 128                         cl = cos(lamtp);
 129                         if (fabs(cl) < tolerance)
 130                             lamtp -= tolerance;
 131                         if( cl < 0 )
 132                             fac = lampp + sin(lampp) * half_pi;
 133                         else
 134                             fac = lampp - sin(lampp) * half_pi;
 135                         for (l = 50; l; --l) {
 136                             lamt = lp_lon + this->m_proj_parm.p22 * sav;
 137                             c = cos(lamt);
 138                             if (fabs(c) < tolerance)
 139                                 lamt -= tolerance;
 140                             xlam = (par.one_es * tanphi * this->m_proj_parm.sa + sin(lamt) * this->m_proj_parm.ca) / c;
 141                             lamdp = atan(xlam) + fac;
 142                             if (fabs(fabs(sav) - fabs(lamdp)) < tolerance)
 143                                 break;
 144                             sav = lamdp;
 145                         }
 146                         if (!l || ++nn >= 3 || (lamdp > this->m_proj_parm.rlm && lamdp < this->m_proj_parm.rlm2))
 147                             break;
 148                         if (lamdp <= this->m_proj_parm.rlm)
 149                             lampp = two_and_half_pi;
 150                         else if (lamdp >= this->m_proj_parm.rlm2)
 151                             lampp = half_pi;
 152                     }
 153                     if (l) {
 154                         sp = sin(lp_lat);
 155                         phidp = aasin((par.one_es * this->m_proj_parm.ca * sp - this->m_proj_parm.sa * cos(lp_lat) *
 156                             sin(lamt)) / sqrt(1. - par.es * sp * sp));
 157                         tanph = log(tan(fourth_pi + .5 * phidp));
 158                         sd = sin(lamdp);
 159                         sdsq = sd * sd;
 160                         s = this->m_proj_parm.p22 * this->m_proj_parm.sa * cos(lamdp) * sqrt((1. + this->m_proj_parm.t * sdsq)
 161                              / ((1. + this->m_proj_parm.w * sdsq) * (1. + this->m_proj_parm.q * sdsq)));
 162                         d = sqrt(this->m_proj_parm.xj * this->m_proj_parm.xj + s * s);
 163                         xy_x = this->m_proj_parm.b * lamdp + this->m_proj_parm.a2 * sin(2. * lamdp) + this->m_proj_parm.a4 *
 164                             sin(lamdp * 4.) - tanph * s / d;
 165                         xy_y = this->m_proj_parm.c1 * sd + this->m_proj_parm.c3 * sin(lamdp * 3.) + tanph * this->m_proj_parm.xj / d;
 166                     } else
 167                         xy_x = xy_y = HUGE_VAL;
 168                 }
 169
 170                 // INVERSE(e_inverse)  ellipsoid
 171                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 172                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
 173                 {
 174                     static const T fourth_pi = detail::fourth_pi<T>();
 175                     static const T half_pi = detail::half_pi<T>();
 176
 177                     int nn;
 178                     T lamt, sdsq, s, lamdp, phidp, sppsq, dd, sd, sl, fac, scl, sav, spp;
 179
 180                     lamdp = xy_x / this->m_proj_parm.b;
 181                     nn = 50;
 182                     do {
 183                         sav = lamdp;
 184                         sd = sin(lamdp);
 185                         sdsq = sd * sd;
 186                         s = this->m_proj_parm.p22 * this->m_proj_parm.sa * cos(lamdp) * sqrt((1. + this->m_proj_parm.t * sdsq)
 187                              / ((1. + this->m_proj_parm.w * sdsq) * (1. + this->m_proj_parm.q * sdsq)));
 188                         lamdp = xy_x + xy_y * s / this->m_proj_parm.xj - this->m_proj_parm.a2 * sin(
 189                             2. * lamdp) - this->m_proj_parm.a4 * sin(lamdp * 4.) - s / this->m_proj_parm.xj * (
 190                             this->m_proj_parm.c1 * sin(lamdp) + this->m_proj_parm.c3 * sin(lamdp * 3.));
 191                         lamdp /= this->m_proj_parm.b;
 192                     } while (fabs(lamdp - sav) >= tolerance && --nn);
 193                     sl = sin(lamdp);
 194                     fac = exp(sqrt(1. + s * s / this->m_proj_parm.xj / this->m_proj_parm.xj) * (xy_y -
 195                         this->m_proj_parm.c1 * sl - this->m_proj_parm.c3 * sin(lamdp * 3.)));
 196                     phidp = 2. * (atan(fac) - fourth_pi);
 197                     dd = sl * sl;
 198                     if (fabs(cos(lamdp)) < tolerance)
 199                         lamdp -= tolerance;
 200                     spp = sin(phidp);
 201                     sppsq = spp * spp;
 202                     lamt = atan(((1. - sppsq * par.rone_es) * tan(lamdp) *
 203                         this->m_proj_parm.ca - spp * this->m_proj_parm.sa * sqrt((1. + this->m_proj_parm.q * dd) * (
 204                         1. - sppsq) - sppsq * this->m_proj_parm.u) / cos(lamdp)) / (1. - sppsq
 205                         * (1. + this->m_proj_parm.u)));
 206                     sl = lamt >= 0. ? 1. : -1.;
 207                     scl = cos(lamdp) >= 0. ? 1. : -1;
 208                     lamt -= half_pi * (1. - scl) * sl;
 209                     lp_lon = lamt - this->m_proj_parm.p22 * lamdp;
 210                     if (fabs(this->m_proj_parm.sa) < tolerance)
 211                         lp_lat = aasin(spp / sqrt(par.one_es * par.one_es + par.es * sppsq));
 212                     else
 213                         lp_lat = atan((tan(lamdp) * cos(lamt) - this->m_proj_parm.ca * sin(lamt)) /
 214                             (par.one_es * this->m_proj_parm.sa));
 215                 }
 216
 217                 static inline std::string get_name()
 218                 {
 219                     return "lsat_ellipsoid";
 220                 }
 221
 222             };
 223
 224             // Space oblique for LANDSAT
 225             template <typename Params, typename Parameters, typename T>
 226             inline void setup_lsat(Params const& params, Parameters& par, par_lsat<T>& proj_parm)
 227             {
 228                 static T const d2r = geometry::math::d2r<T>();
 229                 static T const pi = detail::pi<T>();
 230                 static T const two_pi = detail::two_pi<T>();
 231
 232                 int land, path;
 233                 T lam, alf, esc, ess;
 234
 235                 land = pj_get_param_i<srs::spar::lsat>(params, "lsat", srs::dpar::lsat);
 236                 if (land <= 0 || land > 5)
 237                     BOOST_THROW_EXCEPTION( projection_exception(error_lsat_not_in_range) );
 238
 239                 path = pj_get_param_i<srs::spar::path>(params, "path", srs::dpar::path);
 240                 if (path <= 0 || path > (land <= 3 ? 251 : 233))
 241                     BOOST_THROW_EXCEPTION( projection_exception(error_path_not_in_range) );
 242
 243                 if (land <= 3) {
 244                     par.lam0 = d2r * 128.87 - two_pi / 251. * path;
 245                     proj_parm.p22 = 103.2669323;
 246                     alf = d2r * 99.092;
 247                 } else {
 248                     par.lam0 = d2r * 129.3 - two_pi / 233. * path;
 249                     proj_parm.p22 = 98.8841202;
 250                     alf = d2r * 98.2;
 251                 }
 252                 proj_parm.p22 /= 1440.;
 253                 proj_parm.sa = sin(alf);
 254                 proj_parm.ca = cos(alf);
 255                 if (fabs(proj_parm.ca) < 1e-9)
 256                     proj_parm.ca = 1e-9;
 257                 esc = par.es * proj_parm.ca * proj_parm.ca;
 258                 ess = par.es * proj_parm.sa * proj_parm.sa;
 259                 proj_parm.w = (1. - esc) * par.rone_es;
 260                 proj_parm.w = proj_parm.w * proj_parm.w - 1.;
 261                 proj_parm.q = ess * par.rone_es;
 262                 proj_parm.t = ess * (2. - par.es) * par.rone_es * par.rone_es;
 263                 proj_parm.u = esc * par.rone_es;
 264                 proj_parm.xj = par.one_es * par.one_es * par.one_es;
 265                 proj_parm.rlm = pi * (1. / 248. + .5161290322580645);
 266                 proj_parm.rlm2 = proj_parm.rlm + two_pi;
 267                 proj_parm.a2 = proj_parm.a4 = proj_parm.b = proj_parm.c1 = proj_parm.c3 = 0.;
 268                 seraz0(0., 1., proj_parm);
 269                 for (lam = 9.; lam <= 81.0001; lam += 18.)
 270                     seraz0(lam, 4., proj_parm);
 271                 for (lam = 18; lam <= 72.0001; lam += 18.)
 272                     seraz0(lam, 2., proj_parm);
 273                 seraz0(90., 1., proj_parm);
 274                 proj_parm.a2 /= 30.;
 275                 proj_parm.a4 /= 60.;
 276                 proj_parm.b /= 30.;
 277                 proj_parm.c1 /= 15.;
 278                 proj_parm.c3 /= 45.;
 279             }
 280
 281     }} // namespace detail::lsat
 282     #endif // doxygen
 283
 284     /*!
 285         \brief Space oblique for LANDSAT projection
 286         \ingroup projections
 287         \tparam Geographic latlong point type
 288         \tparam Cartesian xy point type
 289         \tparam Parameters parameter type
 290         \par Projection characteristics
 291          - Cylindrical
 292          - Spheroid
 293          - Ellipsoid
 294         \par Projection parameters
 295          - lsat (integer)
 296          - path (integer)
 297         \par Example
 298         \image html ex_lsat.gif
 299     */
 300     template <typename T, typename Parameters>
 301     struct lsat_ellipsoid : public detail::lsat::base_lsat_ellipsoid<T, Parameters>
 302     {
 303         template <typename Params>
 304         inline lsat_ellipsoid(Params const& params, Parameters & par)
 305         {
 306             detail::lsat::setup_lsat(params, par, this->m_proj_parm);
 307         }
 308     };
 309
 310     #ifndef DOXYGEN_NO_DETAIL
 311     namespace detail
 312     {
 313
 314         // Static projection
 315         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lsat, lsat_ellipsoid)
 316
 317         // Factory entry(s)
 318         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lsat_entry, lsat_ellipsoid)
 319
 320         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(lsat_init)
 321         {
 322             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lsat, lsat_entry)
 323         }
 324
 325     } // namespace detail
 326     #endif // doxygen
 327
 328 } // namespace projections
 329
 330 }} // namespace boost::geometry
 331
 332 #endif // BOOST_GEOMETRY_PROJECTIONS_LSAT_HPP
 333