diff options
Diffstat (limited to 'boost/geometry/srs/projections/proj/omerc.hpp')
-rw-r--r-- | boost/geometry/srs/projections/proj/omerc.hpp | 204 |
1 files changed, 98 insertions, 106 deletions
diff --git a/boost/geometry/srs/projections/proj/omerc.hpp b/boost/geometry/srs/projections/proj/omerc.hpp index 4da6871d13..a4448d182d 100644 --- a/boost/geometry/srs/projections/proj/omerc.hpp +++ b/boost/geometry/srs/projections/proj/omerc.hpp @@ -1,13 +1,9 @@ -#ifndef BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP -#define BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP - -// Boost.Geometry - extensions-gis-projections (based on PROJ4) -// This file is automatically generated. DO NOT EDIT. +// Boost.Geometry - gis-projections (based on PROJ4) // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle. // Use, modification and distribution is subject to the Boost Software License, @@ -19,7 +15,7 @@ // PROJ4 is maintained by Frank Warmerdam // PROJ4 is converted to Boost.Geometry by Barend Gehrels -// Last updated version of proj: 4.9.1 +// Last updated version of proj: 5.0.0 // Original copyright notice: @@ -43,6 +39,9 @@ // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER // DEALINGS IN THE SOFTWARE. +#ifndef BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP +#define BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP + #include <boost/geometry/util/math.hpp> #include <boost/geometry/srs/projections/impl/base_static.hpp> @@ -57,7 +56,7 @@ namespace boost { namespace geometry namespace srs { namespace par4 { - struct omerc {}; + struct omerc {}; // Oblique Mercator }} //namespace srs::par4 @@ -66,9 +65,6 @@ namespace projections #ifndef DOXYGEN_NO_DETAIL namespace detail { namespace omerc { - static const double TOL = 1.e-7; - static const double EPS = 1.e-10; - template <typename T> struct par_omerc { @@ -77,45 +73,45 @@ namespace projections int no_rot; }; + static const double tolerance = 1.e-7; + static const double epsilon = 1.e-10; + // template class, using CRTP to implement forward/inverse - template <typename CalculationType, typename Parameters> - struct base_omerc_ellipsoid : public base_t_fi<base_omerc_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters> + template <typename T, typename Parameters> + struct base_omerc_ellipsoid + : public base_t_fi<base_omerc_ellipsoid<T, Parameters>, T, Parameters> { - - typedef CalculationType geographic_type; - typedef CalculationType cartesian_type; - - par_omerc<CalculationType> m_proj_parm; + par_omerc<T> m_proj_parm; inline base_omerc_ellipsoid(const Parameters& par) - : base_t_fi<base_omerc_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters>(*this, par) {} + : base_t_fi<base_omerc_ellipsoid<T, Parameters>, T, Parameters>(*this, par) + {} // FORWARD(e_forward) ellipsoid // Project coordinates from geographic (lon, lat) to cartesian (x, y) - inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const { - static const CalculationType HALFPI = detail::HALFPI<CalculationType>(); + static const T half_pi = detail::half_pi<T>(); - CalculationType Q, S, T, U, V, temp, u, v; + T s, t, U, V, W, temp, u, v; - if (fabs(fabs(lp_lat) - HALFPI) > EPS) { - Q = this->m_proj_parm.E / pow(pj_tsfn(lp_lat, sin(lp_lat), this->m_par.e), this->m_proj_parm.B); - temp = 1. / Q; - S = .5 * (Q - temp); - T = .5 * (Q + temp); + if (fabs(fabs(lp_lat) - half_pi) > epsilon) { + W = this->m_proj_parm.E / math::pow(pj_tsfn(lp_lat, sin(lp_lat), this->m_par.e), this->m_proj_parm.B); + temp = 1. / W; + s = .5 * (W - temp); + t = .5 * (W + temp); V = sin(this->m_proj_parm.B * lp_lon); - U = (S * this->m_proj_parm.singam - V * this->m_proj_parm.cosgam) / T; - if (fabs(fabs(U) - 1.0) < EPS) - BOOST_THROW_EXCEPTION( projection_exception(-20) ); + U = (s * this->m_proj_parm.singam - V * this->m_proj_parm.cosgam) / t; + if (fabs(fabs(U) - 1.0) < epsilon) { + BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); + } v = 0.5 * this->m_proj_parm.ArB * log((1. - U)/(1. + U)); temp = cos(this->m_proj_parm.B * lp_lon); - if(fabs(temp) < TOL) { - u = this->m_proj_parm.A * lp_lon; - } else { - u = this->m_proj_parm.ArB * atan2((S * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam), temp); - } + if(fabs(temp) < tolerance) { + u = this->m_proj_parm.A * lp_lon; + } else { + u = this->m_proj_parm.ArB * atan2((s * this->m_proj_parm.cosgam + V * this->m_proj_parm.singam), temp); + } } else { v = lp_lat > 0 ? this->m_proj_parm.v_pole_n : this->m_proj_parm.v_pole_s; u = this->m_proj_parm.ArB * lp_lat; @@ -132,11 +128,11 @@ namespace projections // INVERSE(e_inverse) ellipsoid // Project coordinates from cartesian (x, y) to geographic (lon, lat) - inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const { - static const CalculationType HALFPI = detail::HALFPI<CalculationType>(); + static const T half_pi = detail::half_pi<T>(); - CalculationType u, v, Qp, Sp, Tp, Vp, Up; + T u, v, Qp, Sp, Tp, Vp, Up; if (this->m_proj_parm.no_rot) { v = xy_y; @@ -150,13 +146,14 @@ namespace projections Tp = .5 * (Qp + 1. / Qp); Vp = sin(this->m_proj_parm.BrA * u); Up = (Vp * this->m_proj_parm.cosgam + Sp * this->m_proj_parm.singam) / Tp; - if (fabs(fabs(Up) - 1.) < EPS) { + if (fabs(fabs(Up) - 1.) < epsilon) { lp_lon = 0.; - lp_lat = Up < 0. ? -HALFPI : HALFPI; + lp_lat = Up < 0. ? -half_pi : half_pi; } else { lp_lat = this->m_proj_parm.E / sqrt((1. + Up) / (1. - Up)); - if ((lp_lat = pj_phi2(pow(lp_lat, 1. / this->m_proj_parm.B), this->m_par.e)) == HUGE_VAL) - BOOST_THROW_EXCEPTION( projection_exception(-20) ); + if ((lp_lat = pj_phi2(math::pow(lp_lat, T(1) / this->m_proj_parm.B), this->m_par.e)) == HUGE_VAL) { + BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); + } lp_lon = - this->m_proj_parm.rB * atan2((Sp * this->m_proj_parm.cosgam - Vp * this->m_proj_parm.singam), cos(this->m_proj_parm.BrA * u)); } @@ -173,47 +170,46 @@ namespace projections template <typename Parameters, typename T> inline void setup_omerc(Parameters& par, par_omerc<T>& proj_parm) { - static const T FORTPI = detail::FORTPI<T>(); - static const T HALFPI = detail::HALFPI<T>(); - static const T ONEPI = detail::ONEPI<T>(); - static const T TWOPI = detail::TWOPI<T>(); + static const T fourth_pi = detail::fourth_pi<T>(); + static const T half_pi = detail::half_pi<T>(); + static const T pi = detail::pi<T>(); + static const T two_pi = detail::two_pi<T>(); T con, com, cosph0, D, F, H, L, sinph0, p, J, gamma=0, - gamma0, lamc=0, lam1=0, lam2=0, phi1=0, phi2=0, alpha_c=0.0; + gamma0, lamc=0, lam1=0, lam2=0, phi1=0, phi2=0, alpha_c=0; int alp, gam, no_off = 0; - proj_parm.no_rot = pj_param(par.params, "tno_rot").i; - if ((alp = pj_param(par.params, "talpha").i) != 0) - alpha_c = pj_param(par.params, "ralpha").f; - if ((gam = pj_param(par.params, "tgamma").i) != 0) - gamma = pj_param(par.params, "rgamma").f; + proj_parm.no_rot = pj_get_param_b(par.params, "no_rot"); + alp = pj_param_r(par.params, "alpha", alpha_c); + gam = pj_param_r(par.params, "gamma", gamma); if (alp || gam) { - lamc = pj_param(par.params, "rlonc").f; - no_off = - /* For libproj4 compatability */ - pj_param(par.params, "tno_off").i - /* for backward compatibility */ - || pj_param(par.params, "tno_uoff").i; - if( no_off ) - { - /* Mark the parameter as used, so that the pj_get_def() return them */ - pj_param(par.params, "sno_uoff"); - pj_param(par.params, "sno_off"); - } + lamc = pj_get_param_r(par.params, "lonc"); + // NOTE: This is not needed in Boost.Geometry + //no_off = + // /* For libproj4 compatability */ + // pj_param_exists(par.params, "no_off") + // /* for backward compatibility */ + // || pj_param_exists(par.params, "no_uoff"); + //if( no_off ) + //{ + // /* Mark the parameter as used, so that the pj_get_def() return them */ + // pj_get_param_s(par.params, "no_uoff"); + // pj_get_param_s(par.params, "no_off"); + //} } else { - lam1 = pj_param(par.params, "rlon_1").f; - phi1 = pj_param(par.params, "rlat_1").f; - lam2 = pj_param(par.params, "rlon_2").f; - phi2 = pj_param(par.params, "rlat_2").f; - if (fabs(phi1 - phi2) <= TOL || - (con = fabs(phi1)) <= TOL || - fabs(con - HALFPI) <= TOL || - fabs(fabs(par.phi0) - HALFPI) <= TOL || - fabs(fabs(phi2) - HALFPI) <= TOL) - BOOST_THROW_EXCEPTION( projection_exception(-33) ); + lam1 = pj_get_param_r(par.params, "lon_1"); + phi1 = pj_get_param_r(par.params, "lat_1"); + lam2 = pj_get_param_r(par.params, "lon_2"); + phi2 = pj_get_param_r(par.params, "lat_2"); + if (fabs(phi1 - phi2) <= tolerance || + (con = fabs(phi1)) <= tolerance || + fabs(con - half_pi) <= tolerance || + fabs(fabs(par.phi0) - half_pi) <= tolerance || + fabs(fabs(phi2) - half_pi) <= tolerance) + BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_or_alpha_eq_90) ); } com = sqrt(par.one_es); - if (fabs(par.phi0) > EPS) { + if (fabs(par.phi0) > epsilon) { sinph0 = sin(par.phi0); cosph0 = cos(par.phi0); con = 1. - par.es * sinph0 * sinph0; @@ -229,7 +225,7 @@ namespace projections F = -F; } proj_parm.E = F += D; - proj_parm.E *= pow(pj_tsfn(par.phi0, sinph0, par.e), proj_parm.B); + proj_parm.E *= math::pow(pj_tsfn(par.phi0, sinph0, par.e), proj_parm.B); } else { proj_parm.B = 1. / com; proj_parm.A = par.k0; @@ -237,33 +233,29 @@ namespace projections } if (alp || gam) { if (alp) { - gamma0 = asin(sin(alpha_c) / D); + gamma0 = aasin(sin(alpha_c) / D); if (!gam) gamma = alpha_c; } else - alpha_c = asin(D*sin(gamma0 = gamma)); - if ((con = fabs(alpha_c)) <= TOL || - fabs(con - ONEPI) <= TOL || - fabs(fabs(par.phi0) - HALFPI) <= TOL) - BOOST_THROW_EXCEPTION( projection_exception(-32) ); - par.lam0 = lamc - asin(.5 * (F - 1. / F) * + alpha_c = aasin(D*sin(gamma0 = gamma)); + par.lam0 = lamc - aasin(.5 * (F - 1. / F) * tan(gamma0)) / proj_parm.B; } else { - H = pow(pj_tsfn(phi1, sin(phi1), par.e), proj_parm.B); - L = pow(pj_tsfn(phi2, sin(phi2), par.e), proj_parm.B); + H = math::pow(pj_tsfn(phi1, sin(phi1), par.e), proj_parm.B); + L = math::pow(pj_tsfn(phi2, sin(phi2), par.e), proj_parm.B); F = proj_parm.E / H; p = (L - H) / (L + H); J = proj_parm.E * proj_parm.E; J = (J - L * H) / (J + L * H); - if ((con = lam1 - lam2) < -ONEPI) - lam2 -= TWOPI; - else if (con > ONEPI) - lam2 += TWOPI; + if ((con = lam1 - lam2) < -pi) + lam2 -= two_pi; + else if (con > pi) + lam2 += two_pi; par.lam0 = adjlon(.5 * (lam1 + lam2) - atan( J * tan(.5 * proj_parm.B * (lam1 - lam2)) / p) / proj_parm.B); gamma0 = atan(2. * sin(proj_parm.B * adjlon(lam1 - par.lam0)) / (F - 1. / F)); - gamma = alpha_c = asin(D * sin(gamma0)); + gamma = alpha_c = aasin(D * sin(gamma0)); } proj_parm.singam = sin(gamma0); proj_parm.cosgam = cos(gamma0); @@ -274,13 +266,13 @@ namespace projections if (no_off) proj_parm.u_0 = 0; else { - proj_parm.u_0 = fabs(proj_parm.ArB * atan2(sqrt(D * D - 1.), cos(alpha_c))); + proj_parm.u_0 = fabs(proj_parm.ArB * atan(sqrt(D * D - 1.) / cos(alpha_c))); if (par.phi0 < 0.) proj_parm.u_0 = - proj_parm.u_0; } F = 0.5 * gamma0; - proj_parm.v_pole_n = proj_parm.ArB * log(tan(FORTPI - F)); - proj_parm.v_pole_s = proj_parm.ArB * log(tan(FORTPI + F)); + proj_parm.v_pole_n = proj_parm.ArB * log(tan(fourth_pi - F)); + proj_parm.v_pole_s = proj_parm.ArB * log(tan(fourth_pi + F)); } }} // namespace detail::omerc @@ -310,10 +302,10 @@ namespace projections \par Example \image html ex_omerc.gif */ - template <typename CalculationType, typename Parameters> - struct omerc_ellipsoid : public detail::omerc::base_omerc_ellipsoid<CalculationType, Parameters> + template <typename T, typename Parameters> + struct omerc_ellipsoid : public detail::omerc::base_omerc_ellipsoid<T, Parameters> { - inline omerc_ellipsoid(const Parameters& par) : detail::omerc::base_omerc_ellipsoid<CalculationType, Parameters>(par) + inline omerc_ellipsoid(const Parameters& par) : detail::omerc::base_omerc_ellipsoid<T, Parameters>(par) { detail::omerc::setup_omerc(this->m_par, this->m_proj_parm); } @@ -327,20 +319,20 @@ namespace projections BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::omerc, omerc_ellipsoid, omerc_ellipsoid) // Factory entry(s) - template <typename CalculationType, typename Parameters> - class omerc_entry : public detail::factory_entry<CalculationType, Parameters> + template <typename T, typename Parameters> + class omerc_entry : public detail::factory_entry<T, Parameters> { public : - virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + virtual base_v<T, Parameters>* create_new(const Parameters& par) const { - return new base_v_fi<omerc_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + return new base_v_fi<omerc_ellipsoid<T, Parameters>, T, Parameters>(par); } }; - template <typename CalculationType, typename Parameters> - inline void omerc_init(detail::base_factory<CalculationType, Parameters>& factory) + template <typename T, typename Parameters> + inline void omerc_init(detail::base_factory<T, Parameters>& factory) { - factory.add_to_factory("omerc", new omerc_entry<CalculationType, Parameters>); + factory.add_to_factory("omerc", new omerc_entry<T, Parameters>); } } // namespace detail |