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. 6 // Modifications copyright (c) 2017-2018, 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 // 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: 28 29 // The above copyright notice and this permission notice shall be included 30 // in all copies or substantial portions of the Software. 31 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. 39 40 #ifndef BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP 41 #define BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP 42 43 #include <boost/geometry/util/math.hpp> 44 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> 51 52 53 namespace boost { namespace geometry 54 { 55 56 namespace projections 57 { 58 #ifndef DOXYGEN_NO_DETAIL 59 namespace detail { namespace tmerc 60 { 61 62 static const double epsilon10 = 1.e-10; 63 64 template <typename T> FC1()65 inline T FC1() { return 1.; } 66 template <typename T> FC2()67 inline T FC2() { return .5; } 68 template <typename T> FC3()69 inline T FC3() { return .16666666666666666666666666666666666666; } 70 template <typename T> FC4()71 inline T FC4() { return .08333333333333333333333333333333333333; } 72 template <typename T> FC5()73 inline T FC5() { return .05; } 74 template <typename T> FC6()75 inline T FC6() { return .03333333333333333333333333333333333333; } 76 template <typename T> FC7()77 inline T FC7() { return .02380952380952380952380952380952380952; } 78 template <typename T> FC8()79 inline T FC8() { return .01785714285714285714285714285714285714; } 80 81 template <typename T> 82 struct par_tmerc 83 { 84 T esp; 85 T ml0; 86 detail::en<T> en; 87 }; 88 89 // template class, using CRTP to implement forward/inverse 90 template <typename T, typename Parameters> 91 struct base_tmerc_ellipsoid 92 : public base_t_fi<base_tmerc_ellipsoid<T, Parameters>, T, Parameters> 93 { 94 par_tmerc<T> m_proj_parm; 95 base_tmerc_ellipsoidboost::geometry::projections::detail::tmerc::base_tmerc_ellipsoid96 inline base_tmerc_ellipsoid(const Parameters& par) 97 : base_t_fi<base_tmerc_ellipsoid<T, Parameters>, T, Parameters>(*this, par) 98 {} 99 100 // FORWARD(e_forward) ellipse 101 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::tmerc::base_tmerc_ellipsoid102 inline void fwd(T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 103 { 104 static const T half_pi = detail::half_pi<T>(); 105 static const T FC1 = tmerc::FC1<T>(); 106 static const T FC2 = tmerc::FC2<T>(); 107 static const T FC3 = tmerc::FC3<T>(); 108 static const T FC4 = tmerc::FC4<T>(); 109 static const T FC5 = tmerc::FC5<T>(); 110 static const T FC6 = tmerc::FC6<T>(); 111 static const T FC7 = tmerc::FC7<T>(); 112 static const T FC8 = tmerc::FC8<T>(); 113 114 T al, als, n, cosphi, sinphi, t; 115 116 /* 117 * Fail if our longitude is more than 90 degrees from the 118 * central meridian since the results are essentially garbage. 119 * Is error -20 really an appropriate return value? 120 * 121 * http://trac.osgeo.org/proj/ticket/5 122 */ 123 if( lp_lon < -half_pi || lp_lon > half_pi ) 124 { 125 xy_x = HUGE_VAL; 126 xy_y = HUGE_VAL; 127 BOOST_THROW_EXCEPTION( projection_exception(error_lat_or_lon_exceed_limit) ); 128 return; 129 } 130 131 sinphi = sin(lp_lat); 132 cosphi = cos(lp_lat); 133 t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; 134 t *= t; 135 al = cosphi * lp_lon; 136 als = al * al; 137 al /= sqrt(1. - this->m_par.es * sinphi * sinphi); 138 n = this->m_proj_parm.esp * cosphi * cosphi; 139 xy_x = this->m_par.k0 * al * (FC1 + 140 FC3 * als * (1. - t + n + 141 FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t) 142 + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) ) 143 ))); 144 xy_y = this->m_par.k0 * (pj_mlfn(lp_lat, sinphi, cosphi, this->m_proj_parm.en) - this->m_proj_parm.ml0 + 145 sinphi * al * lp_lon * FC2 * ( 1. + 146 FC4 * als * (5. - t + n * (9. + 4. * n) + 147 FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t) 148 + FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) ) 149 )))); 150 } 151 152 // INVERSE(e_inverse) ellipsoid 153 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::tmerc::base_tmerc_ellipsoid154 inline void inv(T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 155 { 156 static const T half_pi = detail::half_pi<T>(); 157 static const T FC1 = tmerc::FC1<T>(); 158 static const T FC2 = tmerc::FC2<T>(); 159 static const T FC3 = tmerc::FC3<T>(); 160 static const T FC4 = tmerc::FC4<T>(); 161 static const T FC5 = tmerc::FC5<T>(); 162 static const T FC6 = tmerc::FC6<T>(); 163 static const T FC7 = tmerc::FC7<T>(); 164 static const T FC8 = tmerc::FC8<T>(); 165 166 T n, con, cosphi, d, ds, sinphi, t; 167 168 lp_lat = pj_inv_mlfn(this->m_proj_parm.ml0 + xy_y / this->m_par.k0, this->m_par.es, this->m_proj_parm.en); 169 if (fabs(lp_lat) >= half_pi) { 170 lp_lat = xy_y < 0. ? -half_pi : half_pi; 171 lp_lon = 0.; 172 } else { 173 sinphi = sin(lp_lat); 174 cosphi = cos(lp_lat); 175 t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; 176 n = this->m_proj_parm.esp * cosphi * cosphi; 177 d = xy_x * sqrt(con = 1. - this->m_par.es * sinphi * sinphi) / this->m_par.k0; 178 con *= t; 179 t *= t; 180 ds = d * d; 181 lp_lat -= (con * ds / (1.-this->m_par.es)) * FC2 * (1. - 182 ds * FC4 * (5. + t * (3. - 9. * n) + n * (1. - 4 * n) - 183 ds * FC6 * (61. + t * (90. - 252. * n + 184 45. * t) + 46. * n 185 - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) ) 186 ))); 187 lp_lon = d*(FC1 - 188 ds*FC3*( 1. + 2.*t + n - 189 ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n 190 - ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) ) 191 ))) / cosphi; 192 } 193 } 194 get_nameboost::geometry::projections::detail::tmerc::base_tmerc_ellipsoid195 static inline std::string get_name() 196 { 197 return "tmerc_ellipsoid"; 198 } 199 200 }; 201 202 // template class, using CRTP to implement forward/inverse 203 template <typename T, typename Parameters> 204 struct base_tmerc_spheroid 205 : public base_t_fi<base_tmerc_spheroid<T, Parameters>, T, Parameters> 206 { 207 par_tmerc<T> m_proj_parm; 208 base_tmerc_spheroidboost::geometry::projections::detail::tmerc::base_tmerc_spheroid209 inline base_tmerc_spheroid(const Parameters& par) 210 : base_t_fi<base_tmerc_spheroid<T, Parameters>, T, Parameters>(*this, par) 211 {} 212 213 // FORWARD(s_forward) sphere 214 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::tmerc::base_tmerc_spheroid215 inline void fwd(T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 216 { 217 static const T half_pi = detail::half_pi<T>(); 218 219 T b, cosphi; 220 221 /* 222 * Fail if our longitude is more than 90 degrees from the 223 * central meridian since the results are essentially garbage. 224 * Is error -20 really an appropriate return value? 225 * 226 * http://trac.osgeo.org/proj/ticket/5 227 */ 228 if( lp_lon < -half_pi || lp_lon > half_pi ) 229 { 230 xy_x = HUGE_VAL; 231 xy_y = HUGE_VAL; 232 BOOST_THROW_EXCEPTION( projection_exception(error_lat_or_lon_exceed_limit) ); 233 return; 234 } 235 236 cosphi = cos(lp_lat); 237 b = cosphi * sin(lp_lon); 238 if (fabs(fabs(b) - 1.) <= epsilon10) 239 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 240 241 xy_x = this->m_proj_parm.ml0 * log((1. + b) / (1. - b)); 242 xy_y = cosphi * cos(lp_lon) / sqrt(1. - b * b); 243 244 b = fabs( xy_y ); 245 if (b >= 1.) { 246 if ((b - 1.) > epsilon10) 247 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 248 else xy_y = 0.; 249 } else 250 xy_y = acos(xy_y); 251 252 if (lp_lat < 0.) 253 xy_y = -xy_y; 254 xy_y = this->m_proj_parm.esp * (xy_y - this->m_par.phi0); 255 } 256 257 // INVERSE(s_inverse) sphere 258 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::tmerc::base_tmerc_spheroid259 inline void inv(T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 260 { 261 T h, g; 262 263 h = exp(xy_x / this->m_proj_parm.esp); 264 g = .5 * (h - 1. / h); 265 h = cos(this->m_par.phi0 + xy_y / this->m_proj_parm.esp); 266 lp_lat = asin(sqrt((1. - h * h) / (1. + g * g))); 267 268 /* Make sure that phi is on the correct hemisphere when false northing is used */ 269 if (xy_y < 0. && -lp_lat+this->m_par.phi0 < 0.0) lp_lat = -lp_lat; 270 271 lp_lon = (g != 0.0 || h != 0.0) ? atan2(g, h) : 0.; 272 } 273 get_nameboost::geometry::projections::detail::tmerc::base_tmerc_spheroid274 static inline std::string get_name() 275 { 276 return "tmerc_spheroid"; 277 } 278 279 }; 280 281 template <typename Parameters, typename T> setup(Parameters & par,par_tmerc<T> & proj_parm)282 inline void setup(Parameters& par, par_tmerc<T>& proj_parm) 283 { 284 if (par.es != 0.0) { 285 proj_parm.en = pj_enfn<T>(par.es); 286 proj_parm.ml0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en); 287 proj_parm.esp = par.es / (1. - par.es); 288 } else { 289 proj_parm.esp = par.k0; 290 proj_parm.ml0 = .5 * proj_parm.esp; 291 } 292 } 293 294 }} // namespace detail::tmerc 295 #endif // doxygen 296 297 /*! 298 \brief Transverse Mercator projection 299 \ingroup projections 300 \tparam Geographic latlong point type 301 \tparam Cartesian xy point type 302 \tparam Parameters parameter type 303 \par Projection characteristics 304 - Cylindrical 305 - Spheroid 306 - Ellipsoid 307 \par Example 308 \image html ex_tmerc.gif 309 */ 310 template <typename T, typename Parameters> 311 struct tmerc_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid<T, Parameters> 312 { 313 template <typename Params> tmerc_ellipsoidboost::geometry::projections::tmerc_ellipsoid314 inline tmerc_ellipsoid(Params const&, Parameters const& par) 315 : detail::tmerc::base_tmerc_ellipsoid<T, Parameters>(par) 316 { 317 detail::tmerc::setup(this->m_par, this->m_proj_parm); 318 } 319 }; 320 321 /*! 322 \brief Transverse Mercator projection 323 \ingroup projections 324 \tparam Geographic latlong point type 325 \tparam Cartesian xy point type 326 \tparam Parameters parameter type 327 \par Projection characteristics 328 - Cylindrical 329 - Spheroid 330 - Ellipsoid 331 \par Example 332 \image html ex_tmerc.gif 333 */ 334 template <typename T, typename Parameters> 335 struct tmerc_spheroid : public detail::tmerc::base_tmerc_spheroid<T, Parameters> 336 { 337 template <typename Params> tmerc_spheroidboost::geometry::projections::tmerc_spheroid338 inline tmerc_spheroid(Params const&, Parameters const& par) 339 : detail::tmerc::base_tmerc_spheroid<T, Parameters>(par) 340 { 341 detail::tmerc::setup(this->m_par, this->m_proj_parm); 342 } 343 }; 344 345 #ifndef DOXYGEN_NO_DETAIL 346 namespace detail 347 { 348 349 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::spar::proj_tmerc,tmerc_spheroid,tmerc_ellipsoid)350 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::spar::proj_tmerc, tmerc_spheroid, tmerc_ellipsoid) 351 352 // Factory entry(s) - dynamic projection 353 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(tmerc_entry, tmerc_spheroid, tmerc_ellipsoid) 354 355 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(tmerc_init) 356 { 357 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(tmerc, tmerc_entry) 358 } 359 360 } // namespace detail 361 #endif // doxygen 362 363 } // namespace projections 364 365 }} // namespace boost::geometry 366 367 #endif // BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP 368 369