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 // 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:
29 // The above copyright notice and this permission notice shall be included
30 // in all copies or substantial portions of the Software.
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.
40 #ifndef BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP
41 #define BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP
43 #include <boost/geometry/util/math.hpp>
45 #include <boost/geometry/srs/projections/impl/base_static.hpp>
46 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
47 #include <boost/geometry/srs/projections/impl/projects.hpp>
48 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
49 #include <boost/geometry/srs/projections/impl/function_overloads.hpp>
50 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
53 namespace boost { namespace geometry
58 #ifndef DOXYGEN_NO_DETAIL
59 namespace detail { namespace tmerc
62 static const double epsilon10 = 1.e-10;
65 inline T FC1() { return 1.; }
67 inline T FC2() { return .5; }
69 inline T FC3() { return .16666666666666666666666666666666666666; }
71 inline T FC4() { return .08333333333333333333333333333333333333; }
73 inline T FC5() { return .05; }
75 inline T FC6() { return .03333333333333333333333333333333333333; }
77 inline T FC7() { return .02380952380952380952380952380952380952; }
79 inline T FC8() { return .01785714285714285714285714285714285714; }
89 template <typename T, typename Parameters>
90 struct base_tmerc_ellipsoid
92 par_tmerc<T> m_proj_parm;
94 // FORWARD(e_forward) ellipse
95 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
96 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
98 static const T half_pi = detail::half_pi<T>();
99 static const T FC1 = tmerc::FC1<T>();
100 static const T FC2 = tmerc::FC2<T>();
101 static const T FC3 = tmerc::FC3<T>();
102 static const T FC4 = tmerc::FC4<T>();
103 static const T FC5 = tmerc::FC5<T>();
104 static const T FC6 = tmerc::FC6<T>();
105 static const T FC7 = tmerc::FC7<T>();
106 static const T FC8 = tmerc::FC8<T>();
108 T al, als, n, cosphi, sinphi, t;
111 * Fail if our longitude is more than 90 degrees from the
112 * central meridian since the results are essentially garbage.
113 * Is error -20 really an appropriate return value?
115 * http://trac.osgeo.org/proj/ticket/5
117 if( lp_lon < -half_pi || lp_lon > half_pi )
121 BOOST_THROW_EXCEPTION( projection_exception(error_lat_or_lon_exceed_limit) );
125 sinphi = sin(lp_lat);
126 cosphi = cos(lp_lat);
127 t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.;
129 al = cosphi * lp_lon;
131 al /= sqrt(1. - par.es * sinphi * sinphi);
132 n = this->m_proj_parm.esp * cosphi * cosphi;
133 xy_x = par.k0 * al * (FC1 +
134 FC3 * als * (1. - t + n +
135 FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t)
136 + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) )
138 xy_y = par.k0 * (pj_mlfn(lp_lat, sinphi, cosphi, this->m_proj_parm.en) - this->m_proj_parm.ml0 +
139 sinphi * al * lp_lon * FC2 * ( 1. +
140 FC4 * als * (5. - t + n * (9. + 4. * n) +
141 FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t)
142 + FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) )
146 // INVERSE(e_inverse) ellipsoid
147 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
148 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
150 static const T half_pi = detail::half_pi<T>();
151 static const T FC1 = tmerc::FC1<T>();
152 static const T FC2 = tmerc::FC2<T>();
153 static const T FC3 = tmerc::FC3<T>();
154 static const T FC4 = tmerc::FC4<T>();
155 static const T FC5 = tmerc::FC5<T>();
156 static const T FC6 = tmerc::FC6<T>();
157 static const T FC7 = tmerc::FC7<T>();
158 static const T FC8 = tmerc::FC8<T>();
160 T n, con, cosphi, d, ds, sinphi, t;
162 lp_lat = pj_inv_mlfn(this->m_proj_parm.ml0 + xy_y / par.k0, par.es, this->m_proj_parm.en);
163 if (fabs(lp_lat) >= half_pi) {
164 lp_lat = xy_y < 0. ? -half_pi : half_pi;
167 sinphi = sin(lp_lat);
168 cosphi = cos(lp_lat);
169 t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.;
170 n = this->m_proj_parm.esp * cosphi * cosphi;
171 d = xy_x * sqrt(con = 1. - par.es * sinphi * sinphi) / par.k0;
175 lp_lat -= (con * ds / (1.-par.es)) * FC2 * (1. -
176 ds * FC4 * (5. + t * (3. - 9. * n) + n * (1. - 4 * n) -
177 ds * FC6 * (61. + t * (90. - 252. * n +
179 - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) )
182 ds*FC3*( 1. + 2.*t + n -
183 ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n
184 - ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) )
189 static inline std::string get_name()
191 return "tmerc_ellipsoid";
196 template <typename T, typename Parameters>
197 struct base_tmerc_spheroid
199 par_tmerc<T> m_proj_parm;
201 // FORWARD(s_forward) sphere
202 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
203 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
205 static const T half_pi = detail::half_pi<T>();
210 * Fail if our longitude is more than 90 degrees from the
211 * central meridian since the results are essentially garbage.
212 * Is error -20 really an appropriate return value?
214 * http://trac.osgeo.org/proj/ticket/5
216 if( lp_lon < -half_pi || lp_lon > half_pi )
220 BOOST_THROW_EXCEPTION( projection_exception(error_lat_or_lon_exceed_limit) );
224 cosphi = cos(lp_lat);
225 b = cosphi * sin(lp_lon);
226 if (fabs(fabs(b) - 1.) <= epsilon10)
227 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
229 xy_x = this->m_proj_parm.ml0 * log((1. + b) / (1. - b));
230 xy_y = cosphi * cos(lp_lon) / sqrt(1. - b * b);
234 if ((b - 1.) > epsilon10)
235 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
242 xy_y = this->m_proj_parm.esp * (xy_y - par.phi0);
245 // INVERSE(s_inverse) sphere
246 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
247 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
251 h = exp(xy_x / this->m_proj_parm.esp);
252 g = .5 * (h - 1. / h);
253 h = cos(par.phi0 + xy_y / this->m_proj_parm.esp);
254 lp_lat = asin(sqrt((1. - h * h) / (1. + g * g)));
256 /* Make sure that phi is on the correct hemisphere when false northing is used */
257 if (xy_y < 0. && -lp_lat+par.phi0 < 0.0) lp_lat = -lp_lat;
259 lp_lon = (g != 0.0 || h != 0.0) ? atan2(g, h) : 0.;
262 static inline std::string get_name()
264 return "tmerc_spheroid";
269 template <typename Parameters, typename T>
270 inline void setup(Parameters const& par, par_tmerc<T>& proj_parm)
273 proj_parm.en = pj_enfn<T>(par.es);
274 proj_parm.ml0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en);
275 proj_parm.esp = par.es / (1. - par.es);
277 proj_parm.esp = par.k0;
278 proj_parm.ml0 = .5 * proj_parm.esp;
282 }} // namespace detail::tmerc
286 \brief Transverse Mercator projection
288 \tparam Geographic latlong point type
289 \tparam Cartesian xy point type
290 \tparam Parameters parameter type
291 \par Projection characteristics
296 \image html ex_tmerc.gif
298 template <typename T, typename Parameters>
299 struct tmerc_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid<T, Parameters>
301 template <typename Params>
302 inline tmerc_ellipsoid(Params const&, Parameters const& par)
304 detail::tmerc::setup(par, this->m_proj_parm);
309 \brief Transverse Mercator projection
311 \tparam Geographic latlong point type
312 \tparam Cartesian xy point type
313 \tparam Parameters parameter type
314 \par Projection characteristics
319 \image html ex_tmerc.gif
321 template <typename T, typename Parameters>
322 struct tmerc_spheroid : public detail::tmerc::base_tmerc_spheroid<T, Parameters>
324 template <typename Params>
325 inline tmerc_spheroid(Params const&, Parameters const& par)
327 detail::tmerc::setup(par, this->m_proj_parm);
331 #ifndef DOXYGEN_NO_DETAIL
336 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_tmerc, tmerc_spheroid, tmerc_ellipsoid)
338 // Factory entry(s) - dynamic projection
339 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(tmerc_entry, tmerc_spheroid, tmerc_ellipsoid)
341 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(tmerc_init)
343 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(tmerc, tmerc_entry)
346 } // namespace detail
349 } // namespace projections
351 }} // namespace boost::geometry
353 #endif // BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP