1 // Boost.Geometry - gis-projections (based on PROJ4)
3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
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.
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)
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
18 // Last updated version of proj: 5.0.0
20 // Original copyright notice:
22 // Copyright (c) 2008 Gerald I. Evenden
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:
31 // The above copyright notice and this permission notice shall be included
32 // in all copies or substantial portions of the Software.
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.
42 /* The code in this file is largly based upon procedures:
44 * Written by: Knud Poder and Karsten Engsager
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.
50 * Modified and used here by permission of Reference Networks
51 * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark
54 #ifndef BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
55 #define BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP
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>
64 #include <boost/math/special_functions/hypot.hpp>
66 namespace boost { namespace geometry
71 #ifndef DOXYGEN_NO_DETAIL
72 namespace detail { namespace etmerc
75 static const int PROJ_ETMERC_ORDER = 6;
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. */
89 inline T log1py(T const& x) { /* Compute log(1+x) accurately */
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;
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;
107 template <typename T>
108 inline T gatg(const T *p1, int len_p1, T const& B) {
110 T h = 0, h1, h2 = 0, cos_2B;
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));
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;
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;
133 for (hi1 = hr1 = hi = 0, hr = *--p; a - p;) {
138 hr = -hr2 + r*hr1 - i*hi1 + *--p;
139 hi = -hi2 + i*hr1 + r*hi1;
141 r = sin_arg_r*cosh_arg_i;
142 i = cos_arg_r*sinh_arg_i;
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;
153 const T* p = a + size;
154 cos_arg_r = cos(arg_r);
158 for (hr1 = 0, hr = *--p; a - p;) {
161 hr = -hr2 + r*hr1 + *--p;
163 return(sin(arg_r)*hr);
166 template <typename T, typename Parameters>
167 struct base_etmerc_ellipsoid
169 par_etmerc<T> m_proj_parm;
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
175 T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
176 T Cn = lp_lat, Ce = lp_lon;
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 */
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));
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);
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 */
197 xy_x = xy_y = HUGE_VAL;
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
204 T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
205 T Cn = xy_y, Ce = xy_x;
208 Cn = (Cn - this->m_proj_parm.Zb)/this->m_proj_parm.Qn;
209 Ce = Ce/this->m_proj_parm.Qn;
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);
215 Ce = atan(sinh(Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - fourth_pi); */
216 /* compl. sph. LAT -> Gaussian LAT, LNG */
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);
228 lp_lat = lp_lon = HUGE_VAL;
231 static inline std::string get_name()
233 return "etmerc_ellipsoid";
238 template <typename Parameters, typename T>
239 inline void setup(Parameters& par, par_etmerc<T>& proj_parm)
244 BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
247 f = par.es / (1 + sqrt(1 - par.es)); /* Replaces: f = 1 - sqrt(1-par.es); */
249 /* third flattening */
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 */
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))))));
262 proj_parm.cgb[1] = np*(7/3.0 + n*( -8/5.0 + n*(-227/45.0 + n*(2704/315.0 +
264 proj_parm.cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0 + n*( 904/315.0 +
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))));
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)));
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));
280 proj_parm.cgb[5] = np*(601676/22275.0 );
281 proj_parm.cbg[5] = np*(444337/155925.0);
283 /* Constants of the projections */
284 /* Transverse Mercator (UTM, ITM, etc) */
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)))));
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))));
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)));
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));
311 proj_parm.utg[5] = np*(-20648693/638668800.0);
312 proj_parm.gtu[5] = np*(212378941/319334400.0);
314 /* Gaussian latitude value of the origin latitude */
315 Z = gatg(proj_parm.cbg, PROJ_ETMERC_ORDER, par.phi0);
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));
322 // Extended Transverse Mercator
323 template <typename Parameters, typename T>
324 inline void setup_etmerc(Parameters& par, par_etmerc<T>& proj_parm)
326 setup(par, proj_parm);
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)
333 static const T pi = detail::pi<T>();
338 BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
341 par.y0 = pj_get_param_b<srs::spar::south>(params, "south", srs::dpar::south) ? 10000000. : 0.;
343 if (pj_param_i<srs::spar::zone>(params, "zone", srs::dpar::zone, zone)) /* zone input ? */
345 if (zone > 0 && zone <= 60)
348 BOOST_THROW_EXCEPTION( projection_exception(error_invalid_utm_zone) );
351 else /* nearest central meridian input */
353 zone = int_floor((adjlon(par.lam0) + pi) * 30. / pi);
359 par.lam0 = (zone + .5) * pi / 30. - pi;
363 setup(par, proj_parm);
366 }} // namespace detail::etmerc
370 \brief Extended Transverse Mercator projection
372 \tparam Geographic latlong point type
373 \tparam Cartesian xy point type
374 \tparam Parameters parameter type
375 \par Projection characteristics
378 \par Projection parameters
379 - lat_ts: Latitude of true scale
380 - lat_0: Latitude of origin
382 \image html ex_etmerc.gif
384 template <typename T, typename Parameters>
385 struct etmerc_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters>
387 template <typename Params>
388 inline etmerc_ellipsoid(Params const& , Parameters & par)
390 detail::etmerc::setup_etmerc(par, this->m_proj_parm);
395 \brief Universal Transverse Mercator (UTM) projection
397 \tparam Geographic latlong point type
398 \tparam Cartesian xy point type
399 \tparam Parameters parameter type
400 \par Projection characteristics
403 \par Projection parameters
404 - zone: UTM Zone (integer)
405 - south: Denotes southern hemisphere UTM zone (boolean)
407 \image html ex_utm.gif
409 template <typename T, typename Parameters>
410 struct utm_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters>
412 template <typename Params>
413 inline utm_ellipsoid(Params const& params, Parameters & par)
415 detail::etmerc::setup_utm(params, par, this->m_proj_parm);
419 #ifndef DOXYGEN_NO_DETAIL
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)
428 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(etmerc_entry, etmerc_ellipsoid)
429 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(utm_entry, utm_ellipsoid)
431 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(etmerc_init)
433 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(etmerc, etmerc_entry);
434 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(utm, utm_entry);
437 } // namespace detail
440 } // namespace projections
442 }} // namespace boost::geometry
444 #endif // BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP