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) 2003, 2006 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 #ifndef BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
43 #define BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
45 #include <boost/geometry/util/math.hpp>
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>
55 namespace boost { namespace geometry
60 #ifndef DOXYGEN_NO_DETAIL
61 namespace detail { namespace omerc
66 T A, B, E, AB, ArB, BrA, rB, singam, cosgam, sinrot, cosrot;
67 T v_pole_n, v_pole_s, u_0;
71 static const double tolerance = 1.e-7;
72 static const double epsilon = 1.e-10;
74 template <typename T, typename Parameters>
75 struct base_omerc_ellipsoid
77 par_omerc<T> m_proj_parm;
79 // FORWARD(e_forward) ellipsoid
80 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
81 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
83 static const T half_pi = detail::half_pi<T>();
85 T s, t, U, V, W, temp, u, v;
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);
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) );
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;
102 u = this->m_proj_parm.ArB * atan2((s * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam), temp);
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;
108 if (this->m_proj_parm.no_rot) {
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;
118 // INVERSE(e_inverse) ellipsoid
119 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
120 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
122 static const T half_pi = detail::half_pi<T>();
124 T u, v, Qp, Sp, Tp, Vp, Up;
126 if (this->m_proj_parm.no_rot) {
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;
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) {
140 lp_lat = Up < 0. ? -half_pi : half_pi;
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) );
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));
151 static inline std::string get_name()
153 return "omerc_ellipsoid";
159 template <typename Params, typename Parameters, typename T>
160 inline void setup_omerc(Params const& params, Parameters & par, par_omerc<T>& proj_parm)
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>();
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;
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);
175 lamc = pj_get_param_r<T, srs::spar::lonc>(params, "lonc", srs::dpar::lonc);
176 // NOTE: This is not needed in Boost.Geometry
178 // /* For libproj4 compatability */
179 // pj_param_exists(par.params, "no_off")
180 // /* for backward compatibility */
181 // || pj_param_exists(par.params, "no_uoff");
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");
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) );
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.)
216 proj_parm.E = F += D;
217 proj_parm.E *= math::pow(pj_tsfn(par.phi0, sinph0, par.e), proj_parm.B);
219 proj_parm.B = 1. / com;
220 proj_parm.A = par.k0;
221 proj_parm.E = D = F = 1.;
225 gamma0 = aasin(sin(alpha_c) / D);
229 alpha_c = aasin(D*sin(gamma0 = gamma));
230 par.lam0 = lamc - aasin(.5 * (F - 1. / F) *
231 tan(gamma0)) / proj_parm.B;
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);
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)
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)) /
247 gamma = alpha_c = aasin(D * sin(gamma0));
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;
258 proj_parm.u_0 = fabs(proj_parm.ArB * atan(sqrt(D * D - 1.) / cos(alpha_c)));
260 proj_parm.u_0 = - proj_parm.u_0;
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));
267 }} // namespace detail::omerc
271 \brief Oblique Mercator projection
273 \tparam Geographic latlong point type
274 \tparam Cartesian xy point type
275 \tparam Parameters parameter type
276 \par Projection characteristics
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)
287 - lat_1: Latitude of first standard parallel (degrees)
289 - lat_2: Latitude of second standard parallel (degrees)
292 \image html ex_omerc.gif
294 template <typename T, typename Parameters>
295 struct omerc_ellipsoid : public detail::omerc::base_omerc_ellipsoid<T, Parameters>
297 template <typename Params>
298 inline omerc_ellipsoid(Params const& params, Parameters & par)
300 detail::omerc::setup_omerc(params, par, this->m_proj_parm);
304 #ifndef DOXYGEN_NO_DETAIL
309 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_omerc, omerc_ellipsoid)
312 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(omerc_entry, omerc_ellipsoid)
314 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(omerc_init)
316 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(omerc, omerc_entry)
319 } // namespace detail
322 } // namespace projections
324 }} // namespace boost::geometry
326 #endif // BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP