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> log1py(T const & x)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> asinhy(T const & x)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> gatg(const T * p1,int len_p1,T const & B)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> clenS(const T * a,int size,T const & arg_r,T const & arg_i,T * R,T * I)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> clens(const T * a,int size,T const & arg_r)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) fwdboost::geometry::projections::detail::etmerc::base_etmerc_ellipsoid173 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) invboost::geometry::projections::detail::etmerc::base_etmerc_ellipsoid202 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 get_nameboost::geometry::projections::detail::etmerc::base_etmerc_ellipsoid231 static inline std::string get_name() 232 { 233 return "etmerc_ellipsoid"; 234 } 235 236 }; 237 238 template <typename Parameters, typename T> setup(Parameters & par,par_etmerc<T> & proj_parm)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> setup_etmerc(Parameters & par,par_etmerc<T> & proj_parm)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> setup_utm(Params const & params,Parameters & par,par_etmerc<T> & proj_parm)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> etmerc_ellipsoidboost::geometry::projections::etmerc_ellipsoid388 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> utm_ellipsoidboost::geometry::projections::utm_ellipsoid413 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 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_etmerc,etmerc_ellipsoid)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 446