boost/geometry/srs/projections/proj/eqdc.hpp

   1 // Boost.Geometry - gis-projections (based on PROJ4)
   2
   3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4
   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.
   8
   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)
  12
  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
  17
  18 // Last updated version of proj: 5.0.0
  19
  20 // Original copyright notice:
  21
  22 // Permission is hereby granted, free of charge, to any person obtaining a
  23 // copy of this software and associated documentation files (the "Software"),
  24 // to deal in the Software without restriction, including without limitation
  25 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  26 // and/or sell copies of the Software, and to permit persons to whom the
  27 // Software is furnished to do so, subject to the following conditions:
  28
  29 // The above copyright notice and this permission notice shall be included
  30 // in all copies or substantial portions of the Software.
  31
  32 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  33 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  34 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  35 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  36 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  37 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  38 // DEALINGS IN THE SOFTWARE.
  39
  40 #ifndef BOOST_GEOMETRY_PROJECTIONS_EQDC_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_EQDC_HPP
  42
  43 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  44 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  45 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  46 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
  47 #include <boost/geometry/srs/projections/impl/pj_msfn.hpp>
  48 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  49 #include <boost/geometry/srs/projections/impl/projects.hpp>
  50
  51 #include <boost/geometry/util/math.hpp>
  52
  53 #include <boost/math/special_functions/hypot.hpp>
  54
  55 namespace boost { namespace geometry
  56 {
  57
  58 namespace projections
  59 {
  60     #ifndef DOXYGEN_NO_DETAIL
  61     namespace detail { namespace eqdc
  62     {
  63
  64             static const double epsilon10 = 1.e-10;
  65
  66             template <typename T>
  67             struct par_eqdc
  68             {
  69                 T    phi1;
  70                 T    phi2;
  71                 T    n;
  72                 T    rho0;
  73                 T    c;
  74                 detail::en<T> en;
  75                 bool ellips;
  76             };
  77
  78             template <typename T, typename Parameters>
  79             struct base_eqdc_ellipsoid
  80             {
  81                 par_eqdc<T> m_proj_parm;
  82
  83                 // FORWARD(e_forward)  sphere & ellipsoid
  84                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  85                 inline void fwd(Parameters const& , T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  86                 {
  87                     T rho = 0.0;
  88
  89                     rho = this->m_proj_parm.c - (this->m_proj_parm.ellips ? pj_mlfn(lp_lat, sin(lp_lat),
  90                         cos(lp_lat), this->m_proj_parm.en) : lp_lat);
  91                     xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
  92                     xy_y = this->m_proj_parm.rho0 - rho * cos(lp_lon);
  93                 }
  94
  95                 // INVERSE(e_inverse)  sphere & ellipsoid
  96                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  97                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
  98                 {
  99                     static T const half_pi = detail::half_pi<T>();
 100
 101                     T rho = 0.0;
 102
 103                     if ((rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho0 - xy_y)) != 0.0 ) {
 104                         if (this->m_proj_parm.n < 0.) {
 105                             rho = -rho;
 106                             xy_x = -xy_x;
 107                             xy_y = -xy_y;
 108                         }
 109                         lp_lat = this->m_proj_parm.c - rho;
 110                         if (this->m_proj_parm.ellips)
 111                             lp_lat = pj_inv_mlfn(lp_lat, par.es, this->m_proj_parm.en);
 112                         lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
 113                     } else {
 114                         lp_lon = 0.;
 115                         lp_lat = this->m_proj_parm.n > 0. ? half_pi : -half_pi;
 116                     }
 117                 }
 118
 119                 static inline std::string get_name()
 120                 {
 121                     return "eqdc_ellipsoid";
 122                 }
 123
 124             };
 125
 126             // Equidistant Conic
 127             template <typename Params, typename Parameters, typename T>
 128             inline void setup_eqdc(Params const& params, Parameters& par, par_eqdc<T>& proj_parm)
 129             {
 130                 T cosphi, sinphi;
 131                 int secant;
 132
 133                 proj_parm.phi1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
 134                 proj_parm.phi2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2);
 135
 136                 if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
 137                     BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
 138
 139                 proj_parm.en = pj_enfn<T>(par.es);
 140
 141                 proj_parm.n = sinphi = sin(proj_parm.phi1);
 142                 cosphi = cos(proj_parm.phi1);
 143                 secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= epsilon10;
 144                 if( (proj_parm.ellips = (par.es > 0.)) ) {
 145                     double ml1, m1;
 146
 147                     m1 = pj_msfn(sinphi, cosphi, par.es);
 148                     ml1 = pj_mlfn(proj_parm.phi1, sinphi, cosphi, proj_parm.en);
 149                     if (secant) { /* secant cone */
 150                         sinphi = sin(proj_parm.phi2);
 151                         cosphi = cos(proj_parm.phi2);
 152                         proj_parm.n = (m1 - pj_msfn(sinphi, cosphi, par.es)) /
 153                             (pj_mlfn(proj_parm.phi2, sinphi, cosphi, proj_parm.en) - ml1);
 154                     }
 155                     proj_parm.c = ml1 + m1 / proj_parm.n;
 156                     proj_parm.rho0 = proj_parm.c - pj_mlfn(par.phi0, sin(par.phi0),
 157                         cos(par.phi0), proj_parm.en);
 158                 } else {
 159                     if (secant)
 160                         proj_parm.n = (cosphi - cos(proj_parm.phi2)) / (proj_parm.phi2 - proj_parm.phi1);
 161                     proj_parm.c = proj_parm.phi1 + cos(proj_parm.phi1) / proj_parm.n;
 162                     proj_parm.rho0 = proj_parm.c - par.phi0;
 163                 }
 164             }
 165
 166     }} // namespace detail::eqdc
 167     #endif // doxygen
 168
 169     /*!
 170         \brief Equidistant Conic projection
 171         \ingroup projections
 172         \tparam Geographic latlong point type
 173         \tparam Cartesian xy point type
 174         \tparam Parameters parameter type
 175         \par Projection characteristics
 176          - Conic
 177          - Spheroid
 178          - Ellipsoid
 179         \par Projection parameters
 180          - lat_1: Latitude of first standard parallel (degrees)
 181          - lat_2: Latitude of second standard parallel (degrees)
 182         \par Example
 183         \image html ex_eqdc.gif
 184     */
 185     template <typename T, typename Parameters>
 186     struct eqdc_ellipsoid : public detail::eqdc::base_eqdc_ellipsoid<T, Parameters>
 187     {
 188         template <typename Params>
 189         inline eqdc_ellipsoid(Params const& params, Parameters const& par)
 190         {
 191             detail::eqdc::setup_eqdc(params, par, this->m_proj_parm);
 192         }
 193     };
 194
 195     #ifndef DOXYGEN_NO_DETAIL
 196     namespace detail
 197     {
 198
 199         // Static projection
 200         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_eqdc, eqdc_ellipsoid)
 201
 202         // Factory entry(s)
 203         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(eqdc_entry, eqdc_ellipsoid)
 204
 205         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(eqdc_init)
 206         {
 207             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(eqdc, eqdc_entry);
 208         }
 209
 210     } // namespace detail
 211     #endif // doxygen
 212
 213 } // namespace projections
 214
 215 }} // namespace boost::geometry
 216
 217 #endif // BOOST_GEOMETRY_PROJECTIONS_EQDC_HPP
 218