Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:35:41 

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 // Copyright (c) 2008   Gerald I. Evenden
0023 
0024 // Permission is hereby granted, free of charge, to any person obtaining a
0025 // copy of this software and associated documentation files (the "Software"),
0026 // to deal in the Software without restriction, including without limitation
0027 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0028 // and/or sell copies of the Software, and to permit persons to whom the
0029 // Software is furnished to do so, subject to the following conditions:
0030 
0031 // The above copyright notice and this permission notice shall be included
0032 // in all copies or substantial portions of the Software.
0033 
0034 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0035 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0036 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0037 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0038 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0039 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0040 // DEALINGS IN THE SOFTWARE.
0041 
0042 /* The code in this file is largly based upon procedures:
0043  *
0044  * Written by: Knud Poder and Karsten Engsager
0045  *
0046  * Based on math from: R.Koenig and K.H. Weise, "Mathematische
0047  * Grundlagen der hoeheren Geodaesie und Kartographie,
0048  * Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951.
0049  *
0050  * Modified and used here by permission of Reference Networks
0051  * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark
0052 */
0053 
0054 #ifndef BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
0055 #define BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
0056 
0057 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0058 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0059 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0060 #include <boost/geometry/srs/projections/impl/function_overloads.hpp>
0061 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0062 #include <boost/geometry/srs/projections/impl/projects.hpp>
0063 
0064 #include <boost/math/special_functions/hypot.hpp>
0065 
0066 namespace boost { namespace geometry
0067 {
0068 
0069 namespace projections
0070 {
0071     #ifndef DOXYGEN_NO_DETAIL
0072     namespace detail { namespace etmerc
0073     {
0074 
0075             static const int PROJ_ETMERC_ORDER = 6;
0076 
0077             template <typename T>
0078             struct par_etmerc
0079             {
0080                 T    Qn;    /* Merid. quad., scaled to the projection */
0081                 T    Zb;    /* Radius vector in polar coord. systems  */
0082                 T    cgb[6]; /* Constants for Gauss -> Geo lat */
0083                 T    cbg[6]; /* Constants for Geo lat -> Gauss */
0084                 T    utg[6]; /* Constants for transv. merc. -> geo */
0085                 T    gtu[6]; /* Constants for geo -> transv. merc. */
0086             };
0087 
0088             template <typename T>
0089             inline T log1py(T const& x) {              /* Compute log(1+x) accurately */
0090                 volatile T
0091                   y = 1 + x,
0092                   z = y - 1;
0093                 /* Here's the explanation for this magic: y = 1 + z, exactly, and z
0094                  * approx x, thus log(y)/z (which is nearly constant near z = 0) returns
0095                  * a good approximation to the true log(1 + x)/x.  The multiplication x *
0096                  * (log(y)/z) introduces little additional error. */
0097                 return z == 0 ? x : x * log(y) / z;
0098             }
0099 
0100             template <typename T>
0101             inline T asinhy(T const& x) {              /* Compute asinh(x) accurately */
0102                 T y = fabs(x);         /* Enforce odd parity */
0103                 y = log1py(y * (1 + y/(boost::math::hypot(1.0, y) + 1)));
0104                 return x < 0 ? -y : y;
0105             }
0106 
0107             template <typename T>
0108             inline T gatg(const T *p1, int len_p1, T const& B) {
0109                 const T *p;
0110                 T h = 0, h1, h2 = 0, cos_2B;
0111 
0112                 cos_2B = 2*cos(2*B);
0113                 for (p = p1 + len_p1, h1 = *--p; p - p1; h2 = h1, h1 = h)
0114                     h = -h2 + cos_2B*h1 + *--p;
0115                 return (B + h*sin(2*B));
0116             }
0117 
0118             /* Complex Clenshaw summation */
0119             template <typename T>
0120             inline T clenS(const T *a, int size, T const& arg_r, T const& arg_i, T *R, T *I) {
0121                 T      r, i, hr, hr1, hr2, hi, hi1, hi2;
0122                 T      sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i;
0123 
0124                 /* arguments */
0125                 const T* p = a + size;
0126                 sin_arg_r  = sin(arg_r);
0127                 cos_arg_r  = cos(arg_r);
0128                 sinh_arg_i = sinh(arg_i);
0129                 cosh_arg_i = cosh(arg_i);
0130                 r          =  2*cos_arg_r*cosh_arg_i;
0131                 i          = -2*sin_arg_r*sinh_arg_i;
0132                 /* summation loop */
0133                 for (hi1 = hr1 = hi = 0, hr = *--p; a - p;) {
0134                     hr2 = hr1;
0135                     hi2 = hi1;
0136                     hr1 = hr;
0137                     hi1 = hi;
0138                     hr  = -hr2 + r*hr1 - i*hi1 + *--p;
0139                     hi  = -hi2 + i*hr1 + r*hi1;
0140                 }
0141                 r   = sin_arg_r*cosh_arg_i;
0142                 i   = cos_arg_r*sinh_arg_i;
0143                 *R  = r*hr - i*hi;
0144                 *I  = r*hi + i*hr;
0145                 return(*R);
0146             }
0147 
0148             /* Real Clenshaw summation */
0149             template <typename T>
0150             inline T clens(const T *a, int size, T const& arg_r) {
0151                 T      r, hr, hr1, hr2, cos_arg_r;
0152 
0153                 const T* p = a + size;
0154                 cos_arg_r  = cos(arg_r);
0155                 r          =  2*cos_arg_r;
0156 
0157                 /* summation loop */
0158                 for (hr1 = 0, hr = *--p; a - p;) {
0159                     hr2 = hr1;
0160                     hr1 = hr;
0161                     hr  = -hr2 + r*hr1 + *--p;
0162                 }
0163                 return(sin(arg_r)*hr);
0164             }
0165 
0166             template <typename T, typename Parameters>
0167             struct base_etmerc_ellipsoid
0168             {
0169                 par_etmerc<T> m_proj_parm;
0170 
0171                 // FORWARD(e_forward)  ellipsoid
0172                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0173                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0174                 {
0175                     T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
0176                     T Cn = lp_lat, Ce = lp_lon;
0177 
0178                     /* ell. LAT, LNG -> Gaussian LAT, LNG */
0179                     Cn  = gatg(this->m_proj_parm.cbg, PROJ_ETMERC_ORDER, Cn);
0180                     /* Gaussian LAT, LNG -> compl. sph. LAT */
0181                     sin_Cn = sin(Cn);
0182                     cos_Cn = cos(Cn);
0183                     sin_Ce = sin(Ce);
0184                     cos_Ce = cos(Ce);
0185 
0186                     Cn     = atan2(sin_Cn, cos_Ce*cos_Cn);
0187                     Ce     = atan2(sin_Ce*cos_Cn, boost::math::hypot(sin_Cn, cos_Cn*cos_Ce));
0188 
0189                     /* compl. sph. N, E -> ell. norm. N, E */
0190                     Ce  = asinhy(tan(Ce));     /* Replaces: Ce  = log(tan(fourth_pi + Ce*0.5)); */
0191                     Cn += clenS(this->m_proj_parm.gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
0192                     Ce += dCe;
0193                     if (fabs(Ce) <= 2.623395162778) {
0194                         xy_y  = this->m_proj_parm.Qn * Cn + this->m_proj_parm.Zb;  /* Northing */
0195                         xy_x  = this->m_proj_parm.Qn * Ce;  /* Easting  */
0196                     } else
0197                         xy_x = xy_y = HUGE_VAL;
0198                 }
0199 
0200                 // INVERSE(e_inverse)  ellipsoid
0201                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0202                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0203                 {
0204                     T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
0205                     T Cn = xy_y, Ce = xy_x;
0206 
0207                     /* normalize N, E */
0208                     Cn = (Cn - this->m_proj_parm.Zb)/this->m_proj_parm.Qn;
0209                     Ce = Ce/this->m_proj_parm.Qn;
0210 
0211                     if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */
0212                         /* norm. N, E -> compl. sph. LAT, LNG */
0213                         Cn += clenS(this->m_proj_parm.utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe);
0214                         Ce += dCe;
0215                         Ce = atan(sinh(Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - fourth_pi); */
0216                         /* compl. sph. LAT -> Gaussian LAT, LNG */
0217                         sin_Cn = sin(Cn);
0218                         cos_Cn = cos(Cn);
0219                         sin_Ce = sin(Ce);
0220                         cos_Ce = cos(Ce);
0221                         Ce     = atan2(sin_Ce, cos_Ce*cos_Cn);
0222                         Cn     = atan2(sin_Cn*cos_Ce, boost::math::hypot(sin_Ce, cos_Ce*cos_Cn));
0223                         /* Gaussian LAT, LNG -> ell. LAT, LNG */
0224                         lp_lat = gatg(this->m_proj_parm.cgb,  PROJ_ETMERC_ORDER, Cn);
0225                         lp_lon = Ce;
0226                     }
0227                     else
0228                         lp_lat = lp_lon = HUGE_VAL;
0229                 }
0230 
0231                 static inline std::string get_name()
0232                 {
0233                     return "etmerc_ellipsoid";
0234                 }
0235 
0236             };
0237 
0238             template <typename Parameters, typename T>
0239             inline void setup(Parameters& par, par_etmerc<T>& proj_parm)
0240             {
0241                 T f, n, np, Z;
0242 
0243                 if (par.es <= 0) {
0244                     BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
0245                 }
0246 
0247                 f = par.es / (1 + sqrt(1 -  par.es)); /* Replaces: f = 1 - sqrt(1-par.es); */
0248 
0249                 /* third flattening */
0250                 np = n = f/(2 - f);
0251 
0252                 /* COEF. OF TRIG SERIES GEO <-> GAUSS */
0253                 /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */
0254                 /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */
0255                 /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */
0256 
0257                 proj_parm.cgb[0] = n*( 2 + n*(-2/3.0  + n*(-2      + n*(116/45.0 + n*(26/45.0 +
0258                             n*(-2854/675.0 ))))));
0259                 proj_parm.cbg[0] = n*(-2 + n*( 2/3.0  + n*( 4/3.0  + n*(-82/45.0 + n*(32/45.0 +
0260                             n*( 4642/4725.0))))));
0261                 np     *= n;
0262                 proj_parm.cgb[1] = np*(7/3.0 + n*( -8/5.0  + n*(-227/45.0 + n*(2704/315.0 +
0263                             n*( 2323/945.0)))));
0264                 proj_parm.cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0  + n*( 904/315.0 +
0265                             n*(-1522/945.0)))));
0266                 np     *= n;
0267                 /* n^5 coeff corrected from 1262/105 -> -1262/105 */
0268                 proj_parm.cgb[2] = np*( 56/15.0  + n*(-136/35.0 + n*(-1262/105.0 +
0269                             n*( 73814/2835.0))));
0270                 proj_parm.cbg[2] = np*(-26/15.0  + n*(  34/21.0 + n*(    8/5.0   +
0271                             n*(-12686/2835.0))));
0272                 np     *= n;
0273                 /* n^5 coeff corrected from 322/35 -> 332/35 */
0274                 proj_parm.cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0)));
0275                 proj_parm.cbg[3] = np*(1237/630.0 + n*( -12/5.0  + n*( -24832/14175.0)));
0276                 np     *= n;
0277                 proj_parm.cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 ));
0278                 proj_parm.cbg[4] = np*(-734/315.0 + n*( 109598/31185.0));
0279                 np     *= n;
0280                 proj_parm.cgb[5] = np*(601676/22275.0 );
0281                 proj_parm.cbg[5] = np*(444337/155925.0);
0282 
0283                 /* Constants of the projections */
0284                 /* Transverse Mercator (UTM, ITM, etc) */
0285                 np = n*n;
0286                 /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */
0287                 proj_parm.Qn = par.k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0)));
0288                 /* coef of trig series */
0289                 /* utg := ell. N, E -> sph. N, E,  KW p194 (65) */
0290                 /* gtu := sph. N, E -> ell. N, E,  KW p196 (69) */
0291                 proj_parm.utg[0] = n*(-0.5  + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 +
0292                             n*(  81/512.0 + n*(-96199/604800.0))))));
0293                 proj_parm.gtu[0] = n*( 0.5  + n*(-2/3.0 + n*(  5/16.0 + n*(41/180.0 +
0294                             n*(-127/288.0 + n*(  7891/37800.0 ))))));
0295                 proj_parm.utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 +
0296                             n*( 1118711/3870720.0)))));
0297                 proj_parm.gtu[1] = np*(13/48.0 + n*(-3/5.0  + n*(557/1440.0 + n*(281/630.0 +
0298                             n*(-1983433/1935360.0)))));
0299                 np      *= n;
0300                 proj_parm.utg[2] = np*(-17/480.0 + n*(  37/840.0 + n*(  209/4480.0  +
0301                             n*( -5569/90720.0 ))));
0302                 proj_parm.gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 +
0303                             n*(167603/181440.0))));
0304                 np      *= n;
0305                 proj_parm.utg[3] = np*(-4397/161280.0 + n*(  11/504.0 + n*( 830251/7257600.0)));
0306                 proj_parm.gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0)));
0307                 np     *= n;
0308                 proj_parm.utg[4] = np*(-4583/161280.0 + n*(  108847/3991680.0));
0309                 proj_parm.gtu[4] = np*(34729/80640.0  + n*(-3418889/1995840.0));
0310                 np     *= n;
0311                 proj_parm.utg[5] = np*(-20648693/638668800.0);
0312                 proj_parm.gtu[5] = np*(212378941/319334400.0);
0313 
0314                 /* Gaussian latitude value of the origin latitude */
0315                 Z = gatg(proj_parm.cbg, PROJ_ETMERC_ORDER, par.phi0);
0316 
0317                 /* Origin northing minus true northing at the origin latitude */
0318                 /* i.e. true northing = N - proj_parm.Zb                         */
0319                 proj_parm.Zb  = - proj_parm.Qn*(Z + clens(proj_parm.gtu, PROJ_ETMERC_ORDER, 2*Z));
0320             }
0321 
0322             // Extended Transverse Mercator
0323             template <typename Parameters, typename T>
0324             inline void setup_etmerc(Parameters& par, par_etmerc<T>& proj_parm)
0325             {
0326                 setup(par, proj_parm);
0327             }
0328 
0329             // Universal Transverse Mercator (UTM)
0330             template <typename Params, typename Parameters, typename T>
0331             inline void setup_utm(Params const& params, Parameters& par, par_etmerc<T>& proj_parm)
0332             {
0333                 static const T pi = detail::pi<T>();
0334 
0335                 int zone;
0336 
0337                 if (par.es == 0.0) {
0338                     BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
0339                 }
0340 
0341                 par.y0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? 10000000. : 0.;
0342                 par.x0 = 500000.;
0343                 if (pj_param_i<srs::spar::zone>(params, "zone", srs::dpar::zone, zone)) /* zone input ? */
0344                 {
0345                     if (zone > 0 && zone <= 60)
0346                         --zone;
0347                     else {
0348                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_utm_zone) );
0349                     }
0350                 }
0351                 else /* nearest central meridian input */
0352                 {
0353                     zone = int_floor((adjlon(par.lam0) + pi) * 30. / pi);
0354                     if (zone < 0)
0355                         zone = 0;
0356                     else if (zone >= 60)
0357                         zone = 59;
0358                 }
0359                 par.lam0 = (zone + .5) * pi / 30. - pi;
0360                 par.k0 = 0.9996;
0361                 par.phi0 = 0.;
0362 
0363                 setup(par, proj_parm);
0364             }
0365 
0366     }} // namespace detail::etmerc
0367     #endif // doxygen
0368 
0369     /*!
0370         \brief Extended Transverse Mercator projection
0371         \ingroup projections
0372         \tparam Geographic latlong point type
0373         \tparam Cartesian xy point type
0374         \tparam Parameters parameter type
0375         \par Projection characteristics
0376          - Cylindrical
0377          - Spheroid
0378         \par Projection parameters
0379          - lat_ts: Latitude of true scale
0380          - lat_0: Latitude of origin
0381         \par Example
0382         \image html ex_etmerc.gif
0383     */
0384     template <typename T, typename Parameters>
0385     struct etmerc_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters>
0386     {
0387         template <typename Params>
0388         inline etmerc_ellipsoid(Params const& , Parameters & par)
0389         {
0390             detail::etmerc::setup_etmerc(par, this->m_proj_parm);
0391         }
0392     };
0393 
0394     /*!
0395         \brief Universal Transverse Mercator (UTM) projection
0396         \ingroup projections
0397         \tparam Geographic latlong point type
0398         \tparam Cartesian xy point type
0399         \tparam Parameters parameter type
0400         \par Projection characteristics
0401          - Cylindrical
0402          - Spheroid
0403         \par Projection parameters
0404          - zone: UTM Zone (integer)
0405          - south: Denotes southern hemisphere UTM zone (boolean)
0406         \par Example
0407         \image html ex_utm.gif
0408     */
0409     template <typename T, typename Parameters>
0410     struct utm_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters>
0411     {
0412         template <typename Params>
0413         inline utm_ellipsoid(Params const& params, Parameters & par)
0414         {
0415             detail::etmerc::setup_utm(params, par, this->m_proj_parm);
0416         }
0417     };
0418 
0419     #ifndef DOXYGEN_NO_DETAIL
0420     namespace detail
0421     {
0422 
0423         // Static projection
0424         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_etmerc, etmerc_ellipsoid)
0425         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_utm, utm_ellipsoid)
0426 
0427         // Factory entry(s)
0428         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(etmerc_entry, etmerc_ellipsoid)
0429         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(utm_entry, utm_ellipsoid)
0430 
0431         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(etmerc_init)
0432         {
0433             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(etmerc, etmerc_entry);
0434             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(utm, utm_entry);
0435         }
0436 
0437     } // namespace detail
0438     #endif // doxygen
0439 
0440 } // namespace projections
0441 
0442 }} // namespace boost::geometry
0443 
0444 #endif // BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
0445
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!