projections/proj/somerc.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.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Permission is hereby granted, free of charge, to any person obtaining a
0023 // copy of this software and associated documentation files (the "Software"),
0024 // to deal in the Software without restriction, including without limitation
0025 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0026 // and/or sell copies of the Software, and to permit persons to whom the
0027 // Software is furnished to do so, subject to the following conditions:
0028
0029 // The above copyright notice and this permission notice shall be included
0030 // in all copies or substantial portions of the Software.
0031
0032 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0033 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0034 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0035 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0036 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0037 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0038 // DEALINGS IN THE SOFTWARE.
0039
0040 #ifndef BOOST_GEOMETRY_PROJECTIONS_SOMERC_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_SOMERC_HPP
0042
0043 #include <boost/geometry/util/math.hpp>
0044
0045 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0046 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0047 #include <boost/geometry/srs/projections/impl/projects.hpp>
0048 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0049 #include <boost/geometry/srs/projections/impl/aasincos.hpp>
0050
0051 namespace boost { namespace geometry
0052 {
0053
0054 namespace projections
0055 {
0056     #ifndef DOXYGEN_NO_DETAIL
0057     namespace detail { namespace somerc
0058     {
0059             static const double epsilon = 1.e-10;
0060             static const int n_iter = 6;
0061
0062             template <typename T>
0063             struct par_somerc
0064             {
0065                 T K, c, hlf_e, kR, cosp0, sinp0;
0066             };
0067
0068             template <typename T, typename Parameters>
0069             struct base_somerc_ellipsoid
0070             {
0071                 par_somerc<T> m_proj_parm;
0072
0073                 // FORWARD(e_forward)
0074                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0075                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0076                 {
0077                     static const T fourth_pi = detail::fourth_pi<T>();
0078                     static const T half_pi = detail::half_pi<T>();
0079
0080                     T phip, lamp, phipp, lampp, sp, cp;
0081
0082                     sp = par.e * sin(lp_lat);
0083                     phip = 2.* atan( exp( this->m_proj_parm.c * (
0084                         log(tan(fourth_pi + 0.5 * lp_lat)) - this->m_proj_parm.hlf_e * log((1. + sp)/(1. - sp)))
0085                         + this->m_proj_parm.K)) - half_pi;
0086                     lamp = this->m_proj_parm.c * lp_lon;
0087                     cp = cos(phip);
0088                     phipp = aasin(this->m_proj_parm.cosp0 * sin(phip) - this->m_proj_parm.sinp0 * cp * cos(lamp));
0089                     lampp = aasin(cp * sin(lamp) / cos(phipp));
0090                     xy_x = this->m_proj_parm.kR * lampp;
0091                     xy_y = this->m_proj_parm.kR * log(tan(fourth_pi + 0.5 * phipp));
0092                 }
0093
0094                 // INVERSE(e_inverse)  ellipsoid & spheroid
0095                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0096                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0097                 {
0098                     static const T fourth_pi = detail::fourth_pi<T>();
0099
0100                     T phip, lamp, phipp, lampp, cp, esp, con, delp;
0101                     int i;
0102
0103                     phipp = 2. * (atan(exp(xy_y / this->m_proj_parm.kR)) - fourth_pi);
0104                     lampp = xy_x / this->m_proj_parm.kR;
0105                     cp = cos(phipp);
0106                     phip = aasin(this->m_proj_parm.cosp0 * sin(phipp) + this->m_proj_parm.sinp0 * cp * cos(lampp));
0107                     lamp = aasin(cp * sin(lampp) / cos(phip));
0108                     con = (this->m_proj_parm.K - log(tan(fourth_pi + 0.5 * phip)))/this->m_proj_parm.c;
0109                     for (i = n_iter; i ; --i) {
0110                         esp = par.e * sin(phip);
0111                         delp = (con + log(tan(fourth_pi + 0.5 * phip)) - this->m_proj_parm.hlf_e *
0112                             log((1. + esp)/(1. - esp)) ) *
0113                             (1. - esp * esp) * cos(phip) * par.rone_es;
0114                         phip -= delp;
0115                         if (fabs(delp) < epsilon)
0116                             break;
0117                     }
0118                     if (i) {
0119                         lp_lat = phip;
0120                         lp_lon = lamp / this->m_proj_parm.c;
0121                     } else {
0122                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0123                     }
0124                 }
0125
0126                 static inline std::string get_name()
0127                 {
0128                     return "somerc_ellipsoid";
0129                 }
0130
0131             };
0132
0133             // Swiss. Obl. Mercator
0134             template <typename Parameters, typename T>
0135             inline void setup_somerc(Parameters const& par, par_somerc<T>& proj_parm)
0136             {
0137                 static const T fourth_pi = detail::fourth_pi<T>();
0138
0139                 T cp, phip0, sp;
0140
0141                 proj_parm.hlf_e = 0.5 * par.e;
0142                 cp = cos(par.phi0);
0143                 cp *= cp;
0144                 proj_parm.c = sqrt(1 + par.es * cp * cp * par.rone_es);
0145                 sp = sin(par.phi0);
0146                 proj_parm.cosp0 = cos( phip0 = aasin(proj_parm.sinp0 = sp / proj_parm.c) );
0147                 sp *= par.e;
0148                 proj_parm.K = log(tan(fourth_pi + 0.5 * phip0)) - proj_parm.c * (
0149                     log(tan(fourth_pi + 0.5 * par.phi0)) - proj_parm.hlf_e *
0150                     log((1. + sp) / (1. - sp)));
0151                 proj_parm.kR = par.k0 * sqrt(par.one_es) / (1. - sp * sp);
0152             }
0153
0154     }} // namespace detail::somerc
0155     #endif // doxygen
0156
0157     /*!
0158         \brief Swiss. Obl. Mercator projection
0159         \ingroup projections
0160         \tparam Geographic latlong point type
0161         \tparam Cartesian xy point type
0162         \tparam Parameters parameter type
0163         \par Projection characteristics
0164          - Cylindrical
0165          - Ellipsoid
0166          - For CH1903
0167         \par Example
0168         \image html ex_somerc.gif
0169     */
0170     template <typename T, typename Parameters>
0171     struct somerc_ellipsoid : public detail::somerc::base_somerc_ellipsoid<T, Parameters>
0172     {
0173         template <typename Params>
0174         inline somerc_ellipsoid(Params const& , Parameters const& par)
0175         {
0176             detail::somerc::setup_somerc(par, this->m_proj_parm);
0177         }
0178     };
0179
0180     #ifndef DOXYGEN_NO_DETAIL
0181     namespace detail
0182     {
0183
0184         // Static projection
0185         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_somerc, somerc_ellipsoid)
0186
0187         // Factory entry(s)
0188         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(somerc_entry, somerc_ellipsoid)
0189
0190         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(somerc_init)
0191         {
0192             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(somerc, somerc_entry)
0193         }
0194
0195     } // namespace detail
0196     #endif // doxygen
0197
0198 } // namespace projections
0199
0200 }} // namespace boost::geometry
0201
0202 #endif // BOOST_GEOMETRY_PROJECTIONS_SOMERC_HPP
0203