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 // 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> 65 inline T FC1() { return 1.; } 66 template <typename T> 67 inline T FC2() { return .5; } 68 template <typename T> 69 inline T FC3() { return .16666666666666666666666666666666666666; } 70 template <typename T> 71 inline T FC4() { return .08333333333333333333333333333333333333; } 72 template <typename T> 73 inline T FC5() { return .05; } 74 template <typename T> 75 inline T FC6() { return .03333333333333333333333333333333333333; } 76 template <typename T> 77 inline T FC7() { return .02380952380952380952380952380952380952; } 78 template <typename T> 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 <typename T, typename Parameters> 90 struct base_tmerc_ellipsoid 91 { 92 par_tmerc<T> m_proj_parm; 93 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 97 { 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>(); 107 108 T al, als, n, cosphi, sinphi, t; 109 110 /* 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? 114 * 115 * http://trac.osgeo.org/proj/ticket/5 116 */ 117 if( lp_lon < -half_pi || lp_lon > half_pi ) 118 { 119 xy_x = HUGE_VAL; 120 xy_y = HUGE_VAL; 121 BOOST_THROW_EXCEPTION( projection_exception(error_lat_or_lon_exceed_limit) ); 122 return; 123 } 124 125 sinphi = sin(lp_lat); 126 cosphi = cos(lp_lat); 127 t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; 128 t *= t; 129 al = cosphi * lp_lon; 130 als = al * al; 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. ) ) 137 ))); 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.) ) 143 )))); 144 } 145 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 149 { 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>(); 159 160 T n, con, cosphi, d, ds, sinphi, t; 161 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; 165 lp_lon = 0.; 166 } else { 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; 172 con *= t; 173 t *= t; 174 ds = d * d; 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 + 178 45. * t) + 46. * n 179 - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) ) 180 ))); 181 lp_lon = d*(FC1 - 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)) ) 185 ))) / cosphi; 186 } 187 } 188 189 static inline std::string get_name() 190 { 191 return "tmerc_ellipsoid"; 192 } 193 194 }; 195 196 template <typename T, typename Parameters> 197 struct base_tmerc_spheroid 198 { 199 par_tmerc<T> m_proj_parm; 200 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 204 { 205 static const T half_pi = detail::half_pi<T>(); 206 207 T b, cosphi; 208 209 /* 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? 213 * 214 * http://trac.osgeo.org/proj/ticket/5 215 */ 216 if( lp_lon < -half_pi || lp_lon > half_pi ) 217 { 218 xy_x = HUGE_VAL; 219 xy_y = HUGE_VAL; 220 BOOST_THROW_EXCEPTION( projection_exception(error_lat_or_lon_exceed_limit) ); 221 return; 222 } 223 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) ); 228 229 xy_x = this->m_proj_parm.ml0 * log((1. + b) / (1. - b)); 230 xy_y = cosphi * cos(lp_lon) / sqrt(1. - b * b); 231 232 b = fabs( xy_y ); 233 if (b >= 1.) { 234 if ((b - 1.) > epsilon10) 235 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 236 else xy_y = 0.; 237 } else 238 xy_y = acos(xy_y); 239 240 if (lp_lat < 0.) 241 xy_y = -xy_y; 242 xy_y = this->m_proj_parm.esp * (xy_y - par.phi0); 243 } 244 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 248 { 249 T h, g; 250 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))); 255 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; 258 259 lp_lon = (g != 0.0 || h != 0.0) ? atan2(g, h) : 0.; 260 } 261 262 static inline std::string get_name() 263 { 264 return "tmerc_spheroid"; 265 } 266 267 }; 268 269 template <typename Parameters, typename T> 270 inline void setup(Parameters const& par, par_tmerc<T>& proj_parm) 271 { 272 if (par.es != 0.0) { 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); 276 } else { 277 proj_parm.esp = par.k0; 278 proj_parm.ml0 = .5 * proj_parm.esp; 279 } 280 } 281 282 }} // namespace detail::tmerc 283 #endif // doxygen 284 285 /*! 286 \brief Transverse Mercator projection 287 \ingroup projections 288 \tparam Geographic latlong point type 289 \tparam Cartesian xy point type 290 \tparam Parameters parameter type 291 \par Projection characteristics 292 - Cylindrical 293 - Spheroid 294 - Ellipsoid 295 \par Example 296 \image html ex_tmerc.gif 297 */ 298 template <typename T, typename Parameters> 299 struct tmerc_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid<T, Parameters> 300 { 301 template <typename Params> 302 inline tmerc_ellipsoid(Params const&, Parameters const& par) 303 { 304 detail::tmerc::setup(par, this->m_proj_parm); 305 } 306 }; 307 308 /*! 309 \brief Transverse Mercator projection 310 \ingroup projections 311 \tparam Geographic latlong point type 312 \tparam Cartesian xy point type 313 \tparam Parameters parameter type 314 \par Projection characteristics 315 - Cylindrical 316 - Spheroid 317 - Ellipsoid 318 \par Example 319 \image html ex_tmerc.gif 320 */ 321 template <typename T, typename Parameters> 322 struct tmerc_spheroid : public detail::tmerc::base_tmerc_spheroid<T, Parameters> 323 { 324 template <typename Params> 325 inline tmerc_spheroid(Params const&, Parameters const& par) 326 { 327 detail::tmerc::setup(par, this->m_proj_parm); 328 } 329 }; 330 331 #ifndef DOXYGEN_NO_DETAIL 332 namespace detail 333 { 334 335 // Static projection 336 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_tmerc, tmerc_spheroid, tmerc_ellipsoid) 337 338 // Factory entry(s) - dynamic projection 339 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(tmerc_entry, tmerc_spheroid, tmerc_ellipsoid) 340 341 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(tmerc_init) 342 { 343 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(tmerc, tmerc_entry) 344 } 345 346 } // namespace detail 347 #endif // doxygen 348 349 } // namespace projections 350 351 }} // namespace boost::geometry 352 353 #endif // BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP 354 355