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) 2003, 2006 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 #ifndef BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP 43 #define BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP 44 45 #include <boost/geometry/util/math.hpp> 46 47 #include <boost/geometry/srs/projections/impl/base_static.hpp> 48 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> 49 #include <boost/geometry/srs/projections/impl/factory_entry.hpp> 50 #include <boost/geometry/srs/projections/impl/pj_param.hpp> 51 #include <boost/geometry/srs/projections/impl/pj_phi2.hpp> 52 #include <boost/geometry/srs/projections/impl/pj_tsfn.hpp> 53 #include <boost/geometry/srs/projections/impl/projects.hpp> 54 55 namespace boost { namespace geometry 56 { 57 58 namespace projections 59 { 60 #ifndef DOXYGEN_NO_DETAIL 61 namespace detail { namespace omerc 62 { 63 template <typename T> 64 struct par_omerc 65 { 66 T A, B, E, AB, ArB, BrA, rB, singam, cosgam, sinrot, cosrot; 67 T v_pole_n, v_pole_s, u_0; 68 bool no_rot; 69 }; 70 71 static const double tolerance = 1.e-7; 72 static const double epsilon = 1.e-10; 73 74 template <typename T, typename Parameters> 75 struct base_omerc_ellipsoid 76 { 77 par_omerc<T> m_proj_parm; 78 79 // FORWARD(e_forward) ellipsoid 80 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::omerc::base_omerc_ellipsoid81 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 82 { 83 static const T half_pi = detail::half_pi<T>(); 84 85 T s, t, U, V, W, temp, u, v; 86 87 if (fabs(fabs(lp_lat) - half_pi) > epsilon) { 88 W = this->m_proj_parm.E / math::pow(pj_tsfn(lp_lat, sin(lp_lat), par.e), this->m_proj_parm.B); 89 temp = 1. / W; 90 s = .5 * (W - temp); 91 t = .5 * (W + temp); 92 V = sin(this->m_proj_parm.B * lp_lon); 93 U = (s * this->m_proj_parm.singam - V * this->m_proj_parm.cosgam) / t; 94 if (fabs(fabs(U) - 1.0) < epsilon) { 95 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 96 } 97 v = 0.5 * this->m_proj_parm.ArB * log((1. - U)/(1. + U)); 98 temp = cos(this->m_proj_parm.B * lp_lon); 99 if(fabs(temp) < tolerance) { 100 u = this->m_proj_parm.A * lp_lon; 101 } else { 102 u = this->m_proj_parm.ArB * atan2((s * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam), temp); 103 } 104 } else { 105 v = lp_lat > 0 ? this->m_proj_parm.v_pole_n : this->m_proj_parm.v_pole_s; 106 u = this->m_proj_parm.ArB * lp_lat; 107 } 108 if (this->m_proj_parm.no_rot) { 109 xy_x = u; 110 xy_y = v; 111 } else { 112 u -= this->m_proj_parm.u_0; 113 xy_x = v * this->m_proj_parm.cosrot + u * this->m_proj_parm.sinrot; 114 xy_y = u * this->m_proj_parm.cosrot - v * this->m_proj_parm.sinrot; 115 } 116 } 117 118 // INVERSE(e_inverse) ellipsoid 119 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::omerc::base_omerc_ellipsoid120 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 121 { 122 static const T half_pi = detail::half_pi<T>(); 123 124 T u, v, Qp, Sp, Tp, Vp, Up; 125 126 if (this->m_proj_parm.no_rot) { 127 v = xy_y; 128 u = xy_x; 129 } else { 130 v = xy_x * this->m_proj_parm.cosrot - xy_y * this->m_proj_parm.sinrot; 131 u = xy_y * this->m_proj_parm.cosrot + xy_x * this->m_proj_parm.sinrot + this->m_proj_parm.u_0; 132 } 133 Qp = exp(- this->m_proj_parm.BrA * v); 134 Sp = .5 * (Qp - 1. / Qp); 135 Tp = .5 * (Qp + 1. / Qp); 136 Vp = sin(this->m_proj_parm.BrA * u); 137 Up = (Vp * this->m_proj_parm.cosgam + Sp * this->m_proj_parm.singam) / Tp; 138 if (fabs(fabs(Up) - 1.) < epsilon) { 139 lp_lon = 0.; 140 lp_lat = Up < 0. ? -half_pi : half_pi; 141 } else { 142 lp_lat = this->m_proj_parm.E / sqrt((1. + Up) / (1. - Up)); 143 if ((lp_lat = pj_phi2(math::pow(lp_lat, T(1) / this->m_proj_parm.B), par.e)) == HUGE_VAL) { 144 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 145 } 146 lp_lon = - this->m_proj_parm.rB * atan2((Sp * this->m_proj_parm.cosgam - 147 Vp * this->m_proj_parm.singam), cos(this->m_proj_parm.BrA * u)); 148 } 149 } 150 get_nameboost::geometry::projections::detail::omerc::base_omerc_ellipsoid151 static inline std::string get_name() 152 { 153 return "omerc_ellipsoid"; 154 } 155 156 }; 157 158 // Oblique Mercator 159 template <typename Params, typename Parameters, typename T> setup_omerc(Params const & params,Parameters & par,par_omerc<T> & proj_parm)160 inline void setup_omerc(Params const& params, Parameters & par, par_omerc<T>& proj_parm) 161 { 162 static const T fourth_pi = detail::fourth_pi<T>(); 163 static const T half_pi = detail::half_pi<T>(); 164 static const T pi = detail::pi<T>(); 165 static const T two_pi = detail::two_pi<T>(); 166 167 T con, com, cosph0, D, F, H, L, sinph0, p, J, gamma=0, 168 gamma0, lamc=0, lam1=0, lam2=0, phi1=0, phi2=0, alpha_c=0; 169 int alp, gam, no_off = 0; 170 171 proj_parm.no_rot = pj_get_param_b<srs::spar::no_rot>(params, "no_rot", srs::dpar::no_rot); 172 alp = pj_param_r<srs::spar::alpha>(params, "alpha", srs::dpar::alpha, alpha_c); 173 gam = pj_param_r<srs::spar::gamma>(params, "gamma", srs::dpar::gamma, gamma); 174 if (alp || gam) { 175 lamc = pj_get_param_r<T, srs::spar::lonc>(params, "lonc", srs::dpar::lonc); 176 // NOTE: This is not needed in Boost.Geometry 177 //no_off = 178 // /* For libproj4 compatability */ 179 // pj_param_exists(par.params, "no_off") 180 // /* for backward compatibility */ 181 // || pj_param_exists(par.params, "no_uoff"); 182 //if( no_off ) 183 //{ 184 // /* Mark the parameter as used, so that the pj_get_def() return them */ 185 // pj_get_param_s(par.params, "no_uoff"); 186 // pj_get_param_s(par.params, "no_off"); 187 //} 188 } else { 189 lam1 = pj_get_param_r<T, srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1); 190 phi1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1); 191 lam2 = pj_get_param_r<T, srs::spar::lon_2>(params, "lon_2", srs::dpar::lon_2); 192 phi2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2); 193 if (fabs(phi1 - phi2) <= tolerance || 194 (con = fabs(phi1)) <= tolerance || 195 fabs(con - half_pi) <= tolerance || 196 fabs(fabs(par.phi0) - half_pi) <= tolerance || 197 fabs(fabs(phi2) - half_pi) <= tolerance) 198 BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_or_alpha_eq_90) ); 199 } 200 com = sqrt(par.one_es); 201 if (fabs(par.phi0) > epsilon) { 202 sinph0 = sin(par.phi0); 203 cosph0 = cos(par.phi0); 204 con = 1. - par.es * sinph0 * sinph0; 205 proj_parm.B = cosph0 * cosph0; 206 proj_parm.B = sqrt(1. + par.es * proj_parm.B * proj_parm.B / par.one_es); 207 proj_parm.A = proj_parm.B * par.k0 * com / con; 208 D = proj_parm.B * com / (cosph0 * sqrt(con)); 209 if ((F = D * D - 1.) <= 0.) 210 F = 0.; 211 else { 212 F = sqrt(F); 213 if (par.phi0 < 0.) 214 F = -F; 215 } 216 proj_parm.E = F += D; 217 proj_parm.E *= math::pow(pj_tsfn(par.phi0, sinph0, par.e), proj_parm.B); 218 } else { 219 proj_parm.B = 1. / com; 220 proj_parm.A = par.k0; 221 proj_parm.E = D = F = 1.; 222 } 223 if (alp || gam) { 224 if (alp) { 225 gamma0 = aasin(sin(alpha_c) / D); 226 if (!gam) 227 gamma = alpha_c; 228 } else 229 alpha_c = aasin(D*sin(gamma0 = gamma)); 230 par.lam0 = lamc - aasin(.5 * (F - 1. / F) * 231 tan(gamma0)) / proj_parm.B; 232 } else { 233 H = math::pow(pj_tsfn(phi1, sin(phi1), par.e), proj_parm.B); 234 L = math::pow(pj_tsfn(phi2, sin(phi2), par.e), proj_parm.B); 235 F = proj_parm.E / H; 236 p = (L - H) / (L + H); 237 J = proj_parm.E * proj_parm.E; 238 J = (J - L * H) / (J + L * H); 239 if ((con = lam1 - lam2) < -pi) 240 lam2 -= two_pi; 241 else if (con > pi) 242 lam2 += two_pi; 243 par.lam0 = adjlon(.5 * (lam1 + lam2) - atan( 244 J * tan(.5 * proj_parm.B * (lam1 - lam2)) / p) / proj_parm.B); 245 gamma0 = atan(2. * sin(proj_parm.B * adjlon(lam1 - par.lam0)) / 246 (F - 1. / F)); 247 gamma = alpha_c = aasin(D * sin(gamma0)); 248 } 249 proj_parm.singam = sin(gamma0); 250 proj_parm.cosgam = cos(gamma0); 251 proj_parm.sinrot = sin(gamma); 252 proj_parm.cosrot = cos(gamma); 253 proj_parm.BrA = 1. / (proj_parm.ArB = proj_parm.A * (proj_parm.rB = 1. / proj_parm.B)); 254 proj_parm.AB = proj_parm.A * proj_parm.B; 255 if (no_off) 256 proj_parm.u_0 = 0; 257 else { 258 proj_parm.u_0 = fabs(proj_parm.ArB * atan(sqrt(D * D - 1.) / cos(alpha_c))); 259 if (par.phi0 < 0.) 260 proj_parm.u_0 = - proj_parm.u_0; 261 } 262 F = 0.5 * gamma0; 263 proj_parm.v_pole_n = proj_parm.ArB * log(tan(fourth_pi - F)); 264 proj_parm.v_pole_s = proj_parm.ArB * log(tan(fourth_pi + F)); 265 } 266 267 }} // namespace detail::omerc 268 #endif // doxygen 269 270 /*! 271 \brief Oblique Mercator projection 272 \ingroup projections 273 \tparam Geographic latlong point type 274 \tparam Cartesian xy point type 275 \tparam Parameters parameter type 276 \par Projection characteristics 277 - Cylindrical 278 - Spheroid 279 - Ellipsoid 280 \par Projection parameters 281 - no_rot: No rotation 282 - alpha: Alpha (degrees) 283 - gamma: Gamma (degrees) 284 - no_off: Only for compatibility with libproj, proj4 (string) 285 - lonc: Longitude (only used if alpha (or gamma) is specified) (degrees) 286 - lon_1 (degrees) 287 - lat_1: Latitude of first standard parallel (degrees) 288 - lon_2 (degrees) 289 - lat_2: Latitude of second standard parallel (degrees) 290 - no_uoff (string) 291 \par Example 292 \image html ex_omerc.gif 293 */ 294 template <typename T, typename Parameters> 295 struct omerc_ellipsoid : public detail::omerc::base_omerc_ellipsoid<T, Parameters> 296 { 297 template <typename Params> omerc_ellipsoidboost::geometry::projections::omerc_ellipsoid298 inline omerc_ellipsoid(Params const& params, Parameters & par) 299 { 300 detail::omerc::setup_omerc(params, par, this->m_proj_parm); 301 } 302 }; 303 304 #ifndef DOXYGEN_NO_DETAIL 305 namespace detail 306 { 307 308 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_omerc,omerc_ellipsoid)309 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_omerc, omerc_ellipsoid) 310 311 // Factory entry(s) 312 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(omerc_entry, omerc_ellipsoid) 313 314 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(omerc_init) 315 { 316 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(omerc, omerc_entry) 317 } 318 319 } // namespace detail 320 #endif // doxygen 321 322 } // namespace projections 323 324 }} // namespace boost::geometry 325 326 #endif // BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP 327 328