projections/proj/omerc.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019, 2022.
0006 // Modifications copyright (c) 2017-2022, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle.
0008 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0009
0010 // Use, modification and distribution is subject to the Boost Software License,
0011 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0012 // http://www.boost.org/LICENSE_1_0.txt)
0013
0014 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0015 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0016 // PROJ4 is maintained by Frank Warmerdam
0017 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0018
0019 // Last updated version of proj: 5.0.0
0020
0021 // Original copyright notice:
0022
0023 // Copyright (c) 2003, 2006   Gerald I. Evenden
0024
0025 // Permission is hereby granted, free of charge, to any person obtaining a
0026 // copy of this software and associated documentation files (the "Software"),
0027 // to deal in the Software without restriction, including without limitation
0028 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0029 // and/or sell copies of the Software, and to permit persons to whom the
0030 // Software is furnished to do so, subject to the following conditions:
0031
0032 // The above copyright notice and this permission notice shall be included
0033 // in all copies or substantial portions of the Software.
0034
0035 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0036 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0037 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0038 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0039 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0040 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0041 // DEALINGS IN THE SOFTWARE.
0042
0043 #ifndef BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
0044 #define BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
0045
0046 #include <boost/geometry/util/math.hpp>
0047
0048 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0049 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0050 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0051 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0052 #include <boost/geometry/srs/projections/impl/pj_phi2.hpp>
0053 #include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
0054 #include <boost/geometry/srs/projections/impl/projects.hpp>
0055
0056 namespace boost { namespace geometry
0057 {
0058
0059 namespace projections
0060 {
0061     #ifndef DOXYGEN_NO_DETAIL
0062     namespace detail { namespace omerc
0063     {
0064             template <typename T>
0065             struct par_omerc
0066             {
0067                 T    A, B, E, AB, ArB, BrA, rB, singam, cosgam, sinrot, cosrot;
0068                 T    v_pole_n, v_pole_s, u_0;
0069                 bool no_rot;
0070             };
0071
0072             static const double tolerance = 1.e-7;
0073             static const double epsilon = 1.e-10;
0074
0075             template <typename T, typename Parameters>
0076             struct base_omerc_ellipsoid
0077             {
0078                 par_omerc<T> m_proj_parm;
0079
0080                 // FORWARD(e_forward)  ellipsoid
0081                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0082                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x,
0083                                 T& xy_y) const
0084                 {
0085                     static const T half_pi = detail::half_pi<T>();
0086
0087                     T  s, t, U, V, W, temp, u, v;
0088
0089                     if (fabs(fabs(lp_lat) - half_pi) > epsilon) {
0090                         W = this->m_proj_parm.E / math::pow(pj_tsfn(lp_lat, sin(lp_lat), par.e),
0091                                                             this->m_proj_parm.B);
0092                         temp = 1. / W;
0093                         s = .5 * (W - temp);
0094                         t = .5 * (W + temp);
0095                         V = sin(this->m_proj_parm.B * lp_lon);
0096                         U = (s * this->m_proj_parm.singam - V * this->m_proj_parm.cosgam) / t;
0097                         if (fabs(fabs(U) - 1.0) < epsilon) {
0098                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0099                         }
0100                         v = 0.5 * this->m_proj_parm.ArB * log((1. - U)/(1. + U));
0101                         temp = cos(this->m_proj_parm.B * lp_lon);
0102                         if(fabs(temp) < tolerance) {
0103                             u = this->m_proj_parm.A * lp_lon;
0104                         } else {
0105                             u = this->m_proj_parm.ArB * atan2((s * this->m_proj_parm.cosgam
0106                                                               + V * this->m_proj_parm.singam), temp);
0107                         }
0108                     } else {
0109                         v = lp_lat > 0 ? this->m_proj_parm.v_pole_n : this->m_proj_parm.v_pole_s;
0110                         u = this->m_proj_parm.ArB * lp_lat;
0111                     }
0112                     if (this->m_proj_parm.no_rot) {
0113                         xy_x = u;
0114                         xy_y = v;
0115                     } else {
0116                         u -= this->m_proj_parm.u_0;
0117                         xy_x = v * this->m_proj_parm.cosrot + u * this->m_proj_parm.sinrot;
0118                         xy_y = u * this->m_proj_parm.cosrot - v * this->m_proj_parm.sinrot;
0119                     }
0120                 }
0121
0122                 // INVERSE(e_inverse)  ellipsoid
0123                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0124                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon,
0125                                 T& lp_lat) const
0126                 {
0127                     static const T half_pi = detail::half_pi<T>();
0128
0129                     T  u, v, Qp, Sp, Tp, Vp, Up;
0130
0131                     if (this->m_proj_parm.no_rot) {
0132                         v = xy_y;
0133                         u = xy_x;
0134                     } else {
0135                         v = xy_x * this->m_proj_parm.cosrot - xy_y * this->m_proj_parm.sinrot;
0136                         u = xy_y * this->m_proj_parm.cosrot + xy_x * this->m_proj_parm.sinrot
0137                             + this->m_proj_parm.u_0;
0138                     }
0139                     Qp = exp(- this->m_proj_parm.BrA * v);
0140                     Sp = .5 * (Qp - 1. / Qp);
0141                     Tp = .5 * (Qp + 1. / Qp);
0142                     Vp = sin(this->m_proj_parm.BrA * u);
0143                     Up = (Vp * this->m_proj_parm.cosgam + Sp * this->m_proj_parm.singam) / Tp;
0144                     if (fabs(fabs(Up) - 1.) < epsilon) {
0145                         lp_lon = 0.;
0146                         lp_lat = Up < 0. ? -half_pi : half_pi;
0147                     } else {
0148                         lp_lat = this->m_proj_parm.E / sqrt((1. + Up) / (1. - Up));
0149                         if ((lp_lat = pj_phi2(math::pow(lp_lat, T(1) / this->m_proj_parm.B), par.e))
0150                             == HUGE_VAL) {
0151                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0152                         }
0153                         lp_lon = - this->m_proj_parm.rB * atan2((Sp * this->m_proj_parm.cosgam -
0154                             Vp * this->m_proj_parm.singam), cos(this->m_proj_parm.BrA * u));
0155                     }
0156                 }
0157
0158                 static inline std::string get_name()
0159                 {
0160                     return "omerc_ellipsoid";
0161                 }
0162
0163             };
0164
0165             // Oblique Mercator
0166             template <typename Params, typename Parameters, typename T>
0167             inline void setup_omerc(Params const& params, Parameters & par, par_omerc<T>& proj_parm)
0168             {
0169                 static const T fourth_pi = detail::fourth_pi<T>();
0170                 static const T half_pi = detail::half_pi<T>();
0171                 static const T pi = detail::pi<T>();
0172                 static const T two_pi = detail::two_pi<T>();
0173
0174                 T con, com, cosph0, D, F, H, L, sinph0, p, J, gamma=0,
0175                   gamma0, lamc=0, lam1=0, lam2=0, phi1=0, phi2=0, alpha_c=0;
0176                 int alp, gam, no_off = 0;
0177
0178                 proj_parm.no_rot = pj_get_param_b<srs::spar::no_rot>(params, "no_rot",
0179                                                                      srs::dpar::no_rot);
0180                 alp = pj_param_r<srs::spar::alpha>(params, "alpha", srs::dpar::alpha, alpha_c);
0181                 gam = pj_param_r<srs::spar::gamma>(params, "gamma", srs::dpar::gamma, gamma);
0182                 if (alp || gam) {
0183                     lamc = pj_get_param_r<T, srs::spar::lonc>(params, "lonc", srs::dpar::lonc);
0184                     // NOTE: This is needed for Hotline Oblique Mercator variant A projection
0185                     no_off = pj_get_param_b<srs::spar::no_off>(params, "no_off", srs::dpar::no_off);
0186                 } else {
0187                     lam1 = pj_get_param_r<T, srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1);
0188                     phi1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
0189                     lam2 = pj_get_param_r<T, srs::spar::lon_2>(params, "lon_2", srs::dpar::lon_2);
0190                     phi2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2);
0191                     if (fabs(phi1 - phi2) <= tolerance ||
0192                         (con = fabs(phi1)) <= tolerance ||
0193                         fabs(con - half_pi) <= tolerance ||
0194                         fabs(fabs(par.phi0) - half_pi) <= tolerance ||
0195                         fabs(fabs(phi2) - half_pi) <= tolerance)
0196                         BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_or_alpha_eq_90) );
0197                 }
0198                 com = sqrt(par.one_es);
0199                 if (fabs(par.phi0) > epsilon) {
0200                     sinph0 = sin(par.phi0);
0201                     cosph0 = cos(par.phi0);
0202                     con = 1. - par.es * sinph0 * sinph0;
0203                     proj_parm.B = cosph0 * cosph0;
0204                     proj_parm.B = sqrt(1. + par.es * proj_parm.B * proj_parm.B / par.one_es);
0205                     proj_parm.A = proj_parm.B * par.k0 * com / con;
0206                     D = proj_parm.B * com / (cosph0 * sqrt(con));
0207                     if ((F = D * D - 1.) <= 0.)
0208                         F = 0.;
0209                     else {
0210                         F = sqrt(F);
0211                         if (par.phi0 < 0.)
0212                             F = -F;
0213                     }
0214                     proj_parm.E = F += D;
0215                     proj_parm.E *= math::pow(pj_tsfn(par.phi0, sinph0, par.e), proj_parm.B);
0216                 } else {
0217                     proj_parm.B = 1. / com;
0218                     proj_parm.A = par.k0;
0219                     proj_parm.E = D = F = 1.;
0220                 }
0221                 if (alp || gam) {
0222                     if (alp) {
0223                         gamma0 = aasin(sin(alpha_c) / D);
0224                         if (!gam)
0225                             gamma = alpha_c;
0226                     } else
0227                         alpha_c = aasin(D*sin(gamma0 = gamma));
0228                     par.lam0 = lamc - aasin(.5 * (F - 1. / F) *
0229                        tan(gamma0)) / proj_parm.B;
0230                 } else {
0231                     H = math::pow(pj_tsfn(phi1, sin(phi1), par.e), proj_parm.B);
0232                     L = math::pow(pj_tsfn(phi2, sin(phi2), par.e), proj_parm.B);
0233                     F = proj_parm.E / H;
0234                     p = (L - H) / (L + H);
0235                     J = proj_parm.E * proj_parm.E;
0236                     J = (J - L * H) / (J + L * H);
0237                     if ((con = lam1 - lam2) < -pi)
0238                         lam2 -= two_pi;
0239                     else if (con > pi)
0240                         lam2 += two_pi;
0241                     par.lam0 = adjlon(.5 * (lam1 + lam2) - atan(
0242                        J * tan(.5 * proj_parm.B * (lam1 - lam2)) / p) / proj_parm.B);
0243                     gamma0 = atan(2. * sin(proj_parm.B * adjlon(lam1 - par.lam0)) /
0244                        (F - 1. / F));
0245                     gamma = alpha_c = aasin(D * sin(gamma0));
0246                 }
0247                 proj_parm.singam = sin(gamma0);
0248                 proj_parm.cosgam = cos(gamma0);
0249                 proj_parm.sinrot = sin(gamma);
0250                 proj_parm.cosrot = cos(gamma);
0251                 proj_parm.BrA = 1. / (proj_parm.ArB = proj_parm.A * (proj_parm.rB = 1. / proj_parm.B));
0252                 proj_parm.AB = proj_parm.A * proj_parm.B;
0253                 if (no_off)
0254                     proj_parm.u_0 = 0;
0255                 else {
0256                     proj_parm.u_0 = fabs(proj_parm.ArB * atan(sqrt(D * D - 1.) / cos(alpha_c)));
0257                     if (par.phi0 < 0.)
0258                         proj_parm.u_0 = - proj_parm.u_0;
0259                 }
0260                 F = 0.5 * gamma0;
0261                 proj_parm.v_pole_n = proj_parm.ArB * log(tan(fourth_pi - F));
0262                 proj_parm.v_pole_s = proj_parm.ArB * log(tan(fourth_pi + F));
0263             }
0264
0265     }} // namespace detail::omerc
0266     #endif // doxygen
0267
0268     /*!
0269         \brief Oblique Mercator projection
0270         \ingroup projections
0271         \tparam Geographic latlong point type
0272         \tparam Cartesian xy point type
0273         \tparam Parameters parameter type
0274         \par Projection characteristics
0275          - Cylindrical
0276          - Spheroid
0277          - Ellipsoid
0278         \par Projection parameters
0279          - no_rot: No rotation
0280          - alpha: Alpha (degrees)
0281          - gamma: Gamma (degrees)
0282          - no_off: Do not offset origin to center of projection
0283             (useful for Hotline Oblique Mercator variant A).
0284          - lonc: Longitude (only used if alpha (or gamma) is specified) (degrees)
0285          - lon_1 (degrees)
0286          - lat_1: Latitude of first standard parallel (degrees)
0287          - lon_2 (degrees)
0288          - lat_2: Latitude of second standard parallel (degrees)
0289          - no_uoff: deprecated (string)
0290         \par Example
0291         \image html ex_omerc.gif
0292     */
0293     template <typename T, typename Parameters>
0294     struct omerc_ellipsoid : public detail::omerc::base_omerc_ellipsoid<T, Parameters>
0295     {
0296         template <typename Params>
0297         inline omerc_ellipsoid(Params const& params, Parameters & par)
0298         {
0299             detail::omerc::setup_omerc(params, par, this->m_proj_parm);
0300         }
0301     };
0302
0303     #ifndef DOXYGEN_NO_DETAIL
0304     namespace detail
0305     {
0306
0307         // Static projection
0308         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_omerc, omerc_ellipsoid)
0309
0310         // Factory entry(s)
0311         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(omerc_entry, omerc_ellipsoid)
0312
0313         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(omerc_init)
0314         {
0315             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(omerc, omerc_entry)
0316         }
0317
0318     } // namespace detail
0319     #endif // doxygen
0320
0321 } // namespace projections
0322
0323 }} // namespace boost::geometry
0324
0325 #endif // BOOST_GEOMETRY_PROJECTIONS_OMERC_HPP
0326