projections/proj/lcc.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_LCC_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_LCC_HPP
0042
0043 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0044 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0045 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0046 #include <boost/geometry/srs/projections/impl/pj_msfn.hpp>
0047 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0048 #include <boost/geometry/srs/projections/impl/pj_phi2.hpp>
0049 #include <boost/geometry/srs/projections/impl/pj_tsfn.hpp>
0050 #include <boost/geometry/srs/projections/impl/projects.hpp>
0051
0052 #include <boost/geometry/util/math.hpp>
0053
0054 #include <boost/math/special_functions/hypot.hpp>
0055
0056 namespace boost { namespace geometry
0057 {
0058
0059 namespace projections
0060 {
0061     #ifndef DOXYGEN_NO_DETAIL
0062     namespace detail { namespace lcc
0063     {
0064             static const double epsilon10 = 1.e-10;
0065
0066             template <typename T>
0067             struct par_lcc
0068             {
0069                 T    phi1;
0070                 T    phi2;
0071                 T    n;
0072                 T    rho0;
0073                 T    c;
0074                 bool ellips;
0075             };
0076
0077             template <typename T, typename Parameters>
0078             struct base_lcc_ellipsoid
0079             {
0080                 par_lcc<T> m_proj_parm;
0081
0082                 // FORWARD(e_forward)  ellipsoid & spheroid
0083                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0084                 inline void fwd(Parameters const& par, T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0085                 {
0086                     static const T fourth_pi = detail::fourth_pi<T>();
0087                     static const T half_pi = detail::half_pi<T>();
0088
0089                     T rho;
0090
0091                     if (fabs(fabs(lp_lat) - half_pi) < epsilon10) {
0092                         if ((lp_lat * this->m_proj_parm.n) <= 0.) {
0093                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0094                         }
0095                         rho = 0.;
0096                     } else {
0097                         rho = this->m_proj_parm.c * (this->m_proj_parm.ellips
0098                             ? math::pow(pj_tsfn(lp_lat, sin(lp_lat), par.e), this->m_proj_parm.n)
0099                             : math::pow(tan(fourth_pi + T(0.5) * lp_lat), -this->m_proj_parm.n));
0100                     }
0101                     lp_lon *= this->m_proj_parm.n;
0102                     xy_x = par.k0 * (rho * sin( lp_lon) );
0103                     xy_y = par.k0 * (this->m_proj_parm.rho0 - rho * cos(lp_lon) );
0104                 }
0105
0106                 // INVERSE(e_inverse)  ellipsoid & spheroid
0107                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0108                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0109                 {
0110                     static const T half_pi = detail::half_pi<T>();
0111
0112                     T rho;
0113
0114                     xy_x /= par.k0;
0115                     xy_y /= par.k0;
0116
0117                     xy_y = this->m_proj_parm.rho0 - xy_y;
0118                     rho = boost::math::hypot(xy_x, xy_y);
0119                     if(rho != 0.0) {
0120                         if (this->m_proj_parm.n < 0.) {
0121                             rho = -rho;
0122                             xy_x = -xy_x;
0123                             xy_y = -xy_y;
0124                         }
0125                         if (this->m_proj_parm.ellips) {
0126                             lp_lat = pj_phi2(math::pow(rho / this->m_proj_parm.c, T(1)/this->m_proj_parm.n), par.e);
0127                             if (lp_lat == HUGE_VAL) {
0128                                 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0129                             }
0130                         } else
0131                             lp_lat = 2. * atan(math::pow(this->m_proj_parm.c / rho, T(1)/this->m_proj_parm.n)) - half_pi;
0132                         lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
0133                     } else {
0134                         lp_lon = 0.;
0135                         lp_lat = this->m_proj_parm.n > 0. ? half_pi : -half_pi;
0136                     }
0137                 }
0138
0139                 static inline std::string get_name()
0140                 {
0141                     return "lcc_ellipsoid";
0142                 }
0143
0144             };
0145
0146             // Lambert Conformal Conic
0147             template <typename Params, typename Parameters, typename T>
0148             inline void setup_lcc(Params const& params, Parameters& par, par_lcc<T>& proj_parm)
0149             {
0150                 static const T fourth_pi = detail::fourth_pi<T>();
0151                 static const T half_pi = detail::half_pi<T>();
0152
0153                 T cosphi, sinphi;
0154                 int secant;
0155
0156                 proj_parm.phi1 = 0.0;
0157                 proj_parm.phi2 = 0.0;
0158                 bool is_phi1_set = pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi1);
0159                 bool is_phi2_set = pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi2);
0160
0161                 // Boost.Geometry specific, set default parameters manually
0162                 if (! is_phi1_set || ! is_phi2_set) {
0163                     bool const use_defaults = ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs);
0164                     if (use_defaults) {
0165                         if (!is_phi1_set) {
0166                             proj_parm.phi1 = 33;
0167                             is_phi1_set = true;
0168                         }
0169                         if (!is_phi2_set) {
0170                             proj_parm.phi2 = 45;
0171                             is_phi2_set = true;
0172                         }
0173                     }
0174                 }
0175
0176                 if (! is_phi2_set) {
0177                     proj_parm.phi2 = proj_parm.phi1;
0178                     if (! pj_param_exists<srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0))
0179                         par.phi0 = proj_parm.phi1;
0180                 }
0181                 if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
0182                     BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
0183
0184                 proj_parm.n = sinphi = sin(proj_parm.phi1);
0185                 cosphi = cos(proj_parm.phi1);
0186                 secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= epsilon10;
0187                 if( (proj_parm.ellips = (par.es != 0.)) ) {
0188                     double ml1, m1;
0189
0190                     par.e = sqrt(par.es); // TODO: Isn't it already set?
0191                     m1 = pj_msfn(sinphi, cosphi, par.es);
0192                     ml1 = pj_tsfn(proj_parm.phi1, sinphi, par.e);
0193                     if (secant) { /* secant cone */
0194                         sinphi = sin(proj_parm.phi2);
0195                         proj_parm.n = log(m1 / pj_msfn(sinphi, cos(proj_parm.phi2), par.es));
0196                         proj_parm.n /= log(ml1 / pj_tsfn(proj_parm.phi2, sinphi, par.e));
0197                     }
0198                     proj_parm.c = (proj_parm.rho0 = m1 * math::pow(ml1, -proj_parm.n) / proj_parm.n);
0199                     proj_parm.rho0 *= (fabs(fabs(par.phi0) - half_pi) < epsilon10) ? T(0) :
0200                         math::pow(pj_tsfn(par.phi0, sin(par.phi0), par.e), proj_parm.n);
0201                 } else {
0202                     if (secant)
0203                         proj_parm.n = log(cosphi / cos(proj_parm.phi2)) /
0204                            log(tan(fourth_pi + .5 * proj_parm.phi2) /
0205                            tan(fourth_pi + .5 * proj_parm.phi1));
0206                     proj_parm.c = cosphi * math::pow(tan(fourth_pi + T(0.5) * proj_parm.phi1), proj_parm.n) / proj_parm.n;
0207                     proj_parm.rho0 = (fabs(fabs(par.phi0) - half_pi) < epsilon10) ? 0. :
0208                         proj_parm.c * math::pow(tan(fourth_pi + T(0.5) * par.phi0), -proj_parm.n);
0209                 }
0210             }
0211
0212     }} // namespace detail::lcc
0213     #endif // doxygen
0214
0215     /*!
0216         \brief Lambert Conformal Conic projection
0217         \ingroup projections
0218         \tparam Geographic latlong point type
0219         \tparam Cartesian xy point type
0220         \tparam Parameters parameter type
0221         \par Projection characteristics
0222          - Conic
0223          - Spheroid
0224          - Ellipsoid
0225         \par Projection parameters
0226          - lat_1: Latitude of first standard parallel (degrees)
0227          - lat_2: Latitude of second standard parallel (degrees)
0228          - lat_0: Latitude of origin
0229         \par Example
0230         \image html ex_lcc.gif
0231     */
0232     template <typename T, typename Parameters>
0233     struct lcc_ellipsoid : public detail::lcc::base_lcc_ellipsoid<T, Parameters>
0234     {
0235         template <typename Params>
0236         inline lcc_ellipsoid(Params const& params, Parameters & par)
0237         {
0238             detail::lcc::setup_lcc(params, par, this->m_proj_parm);
0239         }
0240     };
0241
0242     #ifndef DOXYGEN_NO_DETAIL
0243     namespace detail
0244     {
0245
0246         // Static projection
0247         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_lcc, lcc_ellipsoid)
0248
0249         // Factory entry(s)
0250         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(lcc_entry, lcc_ellipsoid)
0251
0252         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(lcc_init)
0253         {
0254             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(lcc, lcc_entry);
0255         }
0256
0257     } // namespace detail
0258     #endif // doxygen
0259
0260 } // namespace projections
0261
0262 }} // namespace boost::geometry
0263
0264 #endif // BOOST_GEOMETRY_PROJECTIONS_LCC_HPP
0265