boost/geometry/srs/projections/proj/etmerc.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 // Copyright (c) 2008   Gerald I. Evenden
  23
  24 // Permission is hereby granted, free of charge, to any person obtaining a
  25 // copy of this software and associated documentation files (the "Software"),
  26 // to deal in the Software without restriction, including without limitation
  27 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  28 // and/or sell copies of the Software, and to permit persons to whom the
  29 // Software is furnished to do so, subject to the following conditions:
  30
  31 // The above copyright notice and this permission notice shall be included
  32 // in all copies or substantial portions of the Software.
  33
  34 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  35 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  36 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  37 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  38 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  39 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  40 // DEALINGS IN THE SOFTWARE.
  41
  42 /* The code in this file is largly based upon procedures:
  43  *
  44  * Written by: Knud Poder and Karsten Engsager
  45  *
  46  * Based on math from: R.Koenig and K.H. Weise, "Mathematische
  47  * Grundlagen der hoeheren Geodaesie und Kartographie,
  48  * Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951.
  49  *
  50  * Modified and used here by permission of Reference Networks
  51  * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark
  52 */
  53
  54 #ifndef BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
  55 #define BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
  56
  57 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  58 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  59 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  60 #include <boost/geometry/srs/projections/impl/function_overloads.hpp>
  61 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  62 #include <boost/geometry/srs/projections/impl/projects.hpp>
  63
  64 #include <boost/math/special_functions/hypot.hpp>
  65
  66 namespace boost { namespace geometry
  67 {
  68
  69 namespace projections
  70 {
  71     #ifndef DOXYGEN_NO_DETAIL
  72     namespace detail { namespace etmerc
  73     {
  74
  75             static const int PROJ_ETMERC_ORDER = 6;
  76
  77             template <typename T>
  78             struct par_etmerc
  79             {
  80                 T    Qn;    /* Merid. quad., scaled to the projection */
  81                 T    Zb;    /* Radius vector in polar coord. systems  */
  82                 T    cgb[6]; /* Constants for Gauss -> Geo lat */
  83                 T    cbg[6]; /* Constants for Geo lat -> Gauss */
  84                 T    utg[6]; /* Constants for transv. merc. -> geo */
  85                 T    gtu[6]; /* Constants for geo -> transv. merc. */
  86             };
  87
  88             template <typename T>
  89             inline T log1py(T const& x) {              /* Compute log(1+x) accurately */
  90                 volatile T
  91                   y = 1 + x,
  92                   z = y - 1;
  93                 /* Here's the explanation for this magic: y = 1 + z, exactly, and z
  94                  * approx x, thus log(y)/z (which is nearly constant near z = 0) returns
  95                  * a good approximation to the true log(1 + x)/x.  The multiplication x *
  96                  * (log(y)/z) introduces little additional error. */
  97                 return z == 0 ? x : x * log(y) / z;
  98             }
  99
 100             template <typename T>
 101             inline T asinhy(T const& x) {              /* Compute asinh(x) accurately */
 102                 T y = fabs(x);         /* Enforce odd parity */
 103                 y = log1py(y * (1 + y/(boost::math::hypot(1.0, y) + 1)));
 104                 return x < 0 ? -y : y;
 105             }
 106
 107             template <typename T>
 108             inline T gatg(const T *p1, int len_p1, T const& B) {
 109                 const T *p;
 110                 T h = 0, h1, h2 = 0, cos_2B;
 111
 112                 cos_2B = 2*cos(2*B);
 113                 for (p = p1 + len_p1, h1 = *--p; p - p1; h2 = h1, h1 = h)
 114                     h = -h2 + cos_2B*h1 + *--p;
 115                 return (B + h*sin(2*B));
 116             }
 117
 118             /* Complex Clenshaw summation */
 119             template <typename T>
 120             inline T clenS(const T *a, int size, T const& arg_r, T const& arg_i, T *R, T *I) {
 121                 T      r, i, hr, hr1, hr2, hi, hi1, hi2;
 122                 T      sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i;
 123
 124                 /* arguments */
 125                 const T* p = a + size;
 126                 sin_arg_r  = sin(arg_r);
 127                 cos_arg_r  = cos(arg_r);
 128                 sinh_arg_i = sinh(arg_i);
 129                 cosh_arg_i = cosh(arg_i);
 130                 r          =  2*cos_arg_r*cosh_arg_i;
 131                 i          = -2*sin_arg_r*sinh_arg_i;
 132                 /* summation loop */
 133                 for (hi1 = hr1 = hi = 0, hr = *--p; a - p;) {
 134                     hr2 = hr1;
 135                     hi2 = hi1;
 136                     hr1 = hr;
 137                     hi1 = hi;
 138                     hr  = -hr2 + r*hr1 - i*hi1 + *--p;
 139                     hi  = -hi2 + i*hr1 + r*hi1;
 140                 }
 141                 r   = sin_arg_r*cosh_arg_i;
 142                 i   = cos_arg_r*sinh_arg_i;
 143                 *R  = r*hr - i*hi;
 144                 *I  = r*hi + i*hr;
 145                 return(*R);
 146             }
 147
 148             /* Real Clenshaw summation */
 149             template <typename T>
 150             inline T clens(const T *a, int size, T const& arg_r) {
 151                 T      r, hr, hr1, hr2, cos_arg_r;
 152
 153                 const T* p = a + size;
 154                 cos_arg_r  = cos(arg_r);
 155                 r          =  2*cos_arg_r;
 156
 157                 /* summation loop */
 158                 for (hr1 = 0, hr = *--p; a - p;) {
 159                     hr2 = hr1;
 160                     hr1 = hr;
 161                     hr  = -hr2 + r*hr1 + *--p;
 162                 }
 163                 return(sin(arg_r)*hr);
 164             }
 165
 166             template <typename T, typename Parameters>
 167             struct base_etmerc_ellipsoid
 168             {
 169                 par_etmerc<T> m_proj_parm;
 170
 171                 // FORWARD(e_forward)  ellipsoid
 172                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 173                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 174                 {
 175                     T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
 176                     T Cn = lp_lat, Ce = lp_lon;
 177
 178                     /* ell. LAT, LNG -> Gaussian LAT, LNG */
 179                     Cn  = gatg(this->m_proj_parm.cbg, PROJ_ETMERC_ORDER, Cn);
 180                     /* Gaussian LAT, LNG -> compl. sph. LAT */
 181                     sin_Cn = sin(Cn);
 182                     cos_Cn = cos(Cn);
 183                     sin_Ce = sin(Ce);
 184                     cos_Ce = cos(Ce);
 185
 186                     Cn     = atan2(sin_Cn, cos_Ce*cos_Cn);
 187                     Ce     = atan2(sin_Ce*cos_Cn, boost::math::hypot(sin_Cn, cos_Cn*cos_Ce));
 188
 189                     /* compl. sph. N, E -> ell. norm. N, E */
 190                     Ce  = asinhy(tan(Ce));     /* Replaces: Ce  = log(tan(fourth_pi + Ce*0.5)); */
 191                     Cn += clenS(this->m_proj_parm.gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
 192                     Ce += dCe;
 193                     if (fabs(Ce) <= 2.623395162778) {
 194                         xy_y  = this->m_proj_parm.Qn * Cn + this->m_proj_parm.Zb;  /* Northing */
 195                         xy_x  = this->m_proj_parm.Qn * Ce;  /* Easting  */
 196                     } else
 197                         xy_x = xy_y = HUGE_VAL;
 198                 }
 199
 200                 // INVERSE(e_inverse)  ellipsoid
 201                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 202                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
 203                 {
 204                     T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
 205                     T Cn = xy_y, Ce = xy_x;
 206
 207                     /* normalize N, E */
 208                     Cn = (Cn - this->m_proj_parm.Zb)/this->m_proj_parm.Qn;
 209                     Ce = Ce/this->m_proj_parm.Qn;
 210
 211                     if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */
 212                         /* norm. N, E -> compl. sph. LAT, LNG */
 213                         Cn += clenS(this->m_proj_parm.utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
 214                         Ce += dCe;
 215                         Ce = atan(sinh(Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - fourth_pi); */
 216                         /* compl. sph. LAT -> Gaussian LAT, LNG */
 217                         sin_Cn = sin(Cn);
 218                         cos_Cn = cos(Cn);
 219                         sin_Ce = sin(Ce);
 220                         cos_Ce = cos(Ce);
 221                         Ce     = atan2(sin_Ce, cos_Ce*cos_Cn);
 222                         Cn     = atan2(sin_Cn*cos_Ce, boost::math::hypot(sin_Ce, cos_Ce*cos_Cn));
 223                         /* Gaussian LAT, LNG -> ell. LAT, LNG */
 224                         lp_lat = gatg(this->m_proj_parm.cgb,  PROJ_ETMERC_ORDER, Cn);
 225                         lp_lon = Ce;
 226                     }
 227                     else
 228                         lp_lat = lp_lon = HUGE_VAL;
 229                 }
 230
 231                 static inline std::string get_name()
 232                 {
 233                     return "etmerc_ellipsoid";
 234                 }
 235
 236             };
 237
 238             template <typename Parameters, typename T>
 239             inline void setup(Parameters& par, par_etmerc<T>& proj_parm)
 240             {
 241                 T f, n, np, Z;
 242
 243                 if (par.es <= 0) {
 244                     BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
 245                 }
 246
 247                 f = par.es / (1 + sqrt(1 -  par.es)); /* Replaces: f = 1 - sqrt(1-par.es); */
 248
 249                 /* third flattening */
 250                 np = n = f/(2 - f);
 251
 252                 /* COEF. OF TRIG SERIES GEO <-> GAUSS */
 253                 /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */
 254                 /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */
 255                 /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */
 256
 257                 proj_parm.cgb[0] = n*( 2 + n*(-2/3.0  + n*(-2      + n*(116/45.0 + n*(26/45.0 +
 258                             n*(-2854/675.0 ))))));
 259                 proj_parm.cbg[0] = n*(-2 + n*( 2/3.0  + n*( 4/3.0  + n*(-82/45.0 + n*(32/45.0 +
 260                             n*( 4642/4725.0))))));
 261                 np     *= n;
 262                 proj_parm.cgb[1] = np*(7/3.0 + n*( -8/5.0  + n*(-227/45.0 + n*(2704/315.0 +
 263                             n*( 2323/945.0)))));
 264                 proj_parm.cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0  + n*( 904/315.0 +
 265                             n*(-1522/945.0)))));
 266                 np     *= n;
 267                 /* n^5 coeff corrected from 1262/105 -> -1262/105 */
 268                 proj_parm.cgb[2] = np*( 56/15.0  + n*(-136/35.0 + n*(-1262/105.0 +
 269                             n*( 73814/2835.0))));
 270                 proj_parm.cbg[2] = np*(-26/15.0  + n*(  34/21.0 + n*(    8/5.0   +
 271                             n*(-12686/2835.0))));
 272                 np     *= n;
 273                 /* n^5 coeff corrected from 322/35 -> 332/35 */
 274                 proj_parm.cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0)));
 275                 proj_parm.cbg[3] = np*(1237/630.0 + n*( -12/5.0  + n*( -24832/14175.0)));
 276                 np     *= n;
 277                 proj_parm.cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 ));
 278                 proj_parm.cbg[4] = np*(-734/315.0 + n*( 109598/31185.0));
 279                 np     *= n;
 280                 proj_parm.cgb[5] = np*(601676/22275.0 );
 281                 proj_parm.cbg[5] = np*(444337/155925.0);
 282
 283                 /* Constants of the projections */
 284                 /* Transverse Mercator (UTM, ITM, etc) */
 285                 np = n*n;
 286                 /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */
 287                 proj_parm.Qn = par.k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0)));
 288                 /* coef of trig series */
 289                 /* utg := ell. N, E -> sph. N, E,  KW p194 (65) */
 290                 /* gtu := sph. N, E -> ell. N, E,  KW p196 (69) */
 291                 proj_parm.utg[0] = n*(-0.5  + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 +
 292                             n*(  81/512.0 + n*(-96199/604800.0))))));
 293                 proj_parm.gtu[0] = n*( 0.5  + n*(-2/3.0 + n*(  5/16.0 + n*(41/180.0 +
 294                             n*(-127/288.0 + n*(  7891/37800.0 ))))));
 295                 proj_parm.utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 +
 296                             n*( 1118711/3870720.0)))));
 297                 proj_parm.gtu[1] = np*(13/48.0 + n*(-3/5.0  + n*(557/1440.0 + n*(281/630.0 +
 298                             n*(-1983433/1935360.0)))));
 299                 np      *= n;
 300                 proj_parm.utg[2] = np*(-17/480.0 + n*(  37/840.0 + n*(  209/4480.0  +
 301                             n*( -5569/90720.0 ))));
 302                 proj_parm.gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 +
 303                             n*(167603/181440.0))));
 304                 np      *= n;
 305                 proj_parm.utg[3] = np*(-4397/161280.0 + n*(  11/504.0 + n*( 830251/7257600.0)));
 306                 proj_parm.gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0)));
 307                 np     *= n;
 308                 proj_parm.utg[4] = np*(-4583/161280.0 + n*(  108847/3991680.0));
 309                 proj_parm.gtu[4] = np*(34729/80640.0  + n*(-3418889/1995840.0));
 310                 np     *= n;
 311                 proj_parm.utg[5] = np*(-20648693/638668800.0);
 312                 proj_parm.gtu[5] = np*(212378941/319334400.0);
 313
 314                 /* Gaussian latitude value of the origin latitude */
 315                 Z = gatg(proj_parm.cbg, PROJ_ETMERC_ORDER, par.phi0);
 316
 317                 /* Origin northing minus true northing at the origin latitude */
 318                 /* i.e. true northing = N - proj_parm.Zb                         */
 319                 proj_parm.Zb  = - proj_parm.Qn*(Z + clens(proj_parm.gtu, PROJ_ETMERC_ORDER, 2*Z));
 320             }
 321
 322             // Extended Transverse Mercator
 323             template <typename Parameters, typename T>
 324             inline void setup_etmerc(Parameters& par, par_etmerc<T>& proj_parm)
 325             {
 326                 setup(par, proj_parm);
 327             }
 328
 329             // Universal Transverse Mercator (UTM)
 330             template <typename Params, typename Parameters, typename T>
 331             inline void setup_utm(Params const& params, Parameters& par, par_etmerc<T>& proj_parm)
 332             {
 333                 static const T pi = detail::pi<T>();
 334
 335                 int zone;
 336
 337                 if (par.es == 0.0) {
 338                     BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
 339                 }
 340
 341                 par.y0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? 10000000. : 0.;
 342                 par.x0 = 500000.;
 343                 if (pj_param_i<srs::spar::zone>(params, "zone", srs::dpar::zone, zone)) /* zone input ? */
 344                 {
 345                     if (zone > 0 && zone <= 60)
 346                         --zone;
 347                     else {
 348                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_utm_zone) );
 349                     }
 350                 }
 351                 else /* nearest central meridian input */
 352                 {
 353                     zone = int_floor((adjlon(par.lam0) + pi) * 30. / pi);
 354                     if (zone < 0)
 355                         zone = 0;
 356                     else if (zone >= 60)
 357                         zone = 59;
 358                 }
 359                 par.lam0 = (zone + .5) * pi / 30. - pi;
 360                 par.k0 = 0.9996;
 361                 par.phi0 = 0.;
 362
 363                 setup(par, proj_parm);
 364             }
 365
 366     }} // namespace detail::etmerc
 367     #endif // doxygen
 368
 369     /*!
 370         \brief Extended Transverse Mercator projection
 371         \ingroup projections
 372         \tparam Geographic latlong point type
 373         \tparam Cartesian xy point type
 374         \tparam Parameters parameter type
 375         \par Projection characteristics
 376          - Cylindrical
 377          - Spheroid
 378         \par Projection parameters
 379          - lat_ts: Latitude of true scale
 380          - lat_0: Latitude of origin
 381         \par Example
 382         \image html ex_etmerc.gif
 383     */
 384     template <typename T, typename Parameters>
 385     struct etmerc_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters>
 386     {
 387         template <typename Params>
 388         inline etmerc_ellipsoid(Params const& , Parameters & par)
 389         {
 390             detail::etmerc::setup_etmerc(par, this->m_proj_parm);
 391         }
 392     };
 393
 394     /*!
 395         \brief Universal Transverse Mercator (UTM) projection
 396         \ingroup projections
 397         \tparam Geographic latlong point type
 398         \tparam Cartesian xy point type
 399         \tparam Parameters parameter type
 400         \par Projection characteristics
 401          - Cylindrical
 402          - Spheroid
 403         \par Projection parameters
 404          - zone: UTM Zone (integer)
 405          - south: Denotes southern hemisphere UTM zone (boolean)
 406         \par Example
 407         \image html ex_utm.gif
 408     */
 409     template <typename T, typename Parameters>
 410     struct utm_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters>
 411     {
 412         template <typename Params>
 413         inline utm_ellipsoid(Params const& params, Parameters & par)
 414         {
 415             detail::etmerc::setup_utm(params, par, this->m_proj_parm);
 416         }
 417     };
 418
 419     #ifndef DOXYGEN_NO_DETAIL
 420     namespace detail
 421     {
 422
 423         // Static projection
 424         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_etmerc, etmerc_ellipsoid)
 425         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_utm, utm_ellipsoid)
 426
 427         // Factory entry(s)
 428         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(etmerc_entry, etmerc_ellipsoid)
 429         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(utm_entry, utm_ellipsoid)
 430
 431         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(etmerc_init)
 432         {
 433             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(etmerc, etmerc_entry);
 434             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(utm, utm_entry);
 435         }
 436
 437     } // namespace detail
 438     #endif // doxygen
 439
 440 } // namespace projections
 441
 442 }} // namespace boost::geometry
 443
 444 #endif // BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
 445