boost/geometry/srs/projections/proj/ortho.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_ORTHO_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_ORTHO_HPP
  42
  43 #include <boost/config.hpp>
  44 #include <boost/geometry/util/math.hpp>
  45 #include <boost/math/special_functions/hypot.hpp>
  46
  47 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  48 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  49 #include <boost/geometry/srs/projections/impl/projects.hpp>
  50 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  51
  52 namespace boost { namespace geometry
  53 {
  54
  55 namespace projections
  56 {
  57     #ifndef DOXYGEN_NO_DETAIL
  58     namespace detail { namespace ortho
  59     {
  60
  61             enum mode_type {
  62                 n_pole = 0,
  63                 s_pole = 1,
  64                 equit  = 2,
  65                 obliq  = 3
  66             };
  67
  68             template <typename T>
  69             struct par_ortho
  70             {
  71                 T   sinph0;
  72                 T   cosph0;
  73                 mode_type mode;
  74             };
  75
  76             static const double epsilon10 = 1.e-10;
  77
  78             template <typename T, typename Parameters>
  79             struct base_ortho_spheroid
  80             {
  81                 par_ortho<T> m_proj_parm;
  82
  83                 // FORWARD(s_forward)  spheroid
  84                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  85                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  86                 {
  87                     static const T half_pi = detail::half_pi<T>();
  88
  89                     T coslam, cosphi, sinphi;
  90
  91                     cosphi = cos(lp_lat);
  92                     coslam = cos(lp_lon);
  93                     switch (this->m_proj_parm.mode) {
  94                     case equit:
  95                         if (cosphi * coslam < - epsilon10) {
  96                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
  97                         }
  98                         xy_y = sin(lp_lat);
  99                         break;
 100                     case obliq:
 101                         if (this->m_proj_parm.sinph0 * (sinphi = sin(lp_lat)) +
 102                            this->m_proj_parm.cosph0 * cosphi * coslam < - epsilon10) {
 103                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 104                         }
 105                         xy_y = this->m_proj_parm.cosph0 * sinphi - this->m_proj_parm.sinph0 * cosphi * coslam;
 106                         break;
 107                     case n_pole:
 108                         coslam = - coslam;
 109                         BOOST_FALLTHROUGH;
 110                     case s_pole:
 111                         if (fabs(lp_lat - par.phi0) - epsilon10 > half_pi) {
 112                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 113                         }
 114                         xy_y = cosphi * coslam;
 115                         break;
 116                     }
 117                     xy_x = cosphi * sin(lp_lon);
 118                 }
 119
 120                 // INVERSE(s_inverse)  spheroid
 121                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 122                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 123                 {
 124                     static const T half_pi = detail::half_pi<T>();
 125
 126                     T rh, cosc, sinc;
 127
 128                     if ((sinc = (rh = boost::math::hypot(xy_x, xy_y))) > 1.) {
 129                         if ((sinc - 1.) > epsilon10) {
 130                             BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
 131                         }
 132                         sinc = 1.;
 133                     }
 134                     cosc = sqrt(1. - sinc * sinc); /* in this range OK */
 135                     if (fabs(rh) <= epsilon10) {
 136                         lp_lat = par.phi0;
 137                         lp_lon = 0.0;
 138                     } else {
 139                         switch (this->m_proj_parm.mode) {
 140                         case n_pole:
 141                             xy_y = -xy_y;
 142                             lp_lat = acos(sinc);
 143                             break;
 144                         case s_pole:
 145                             lp_lat = - acos(sinc);
 146                             break;
 147                         case equit:
 148                             lp_lat = xy_y * sinc / rh;
 149                             xy_x *= sinc;
 150                             xy_y = cosc * rh;
 151                             goto sinchk;
 152                         case obliq:
 153                             lp_lat = cosc * this->m_proj_parm.sinph0 + xy_y * sinc * this->m_proj_parm.cosph0 /rh;
 154                             xy_y = (cosc - this->m_proj_parm.sinph0 * lp_lat) * rh;
 155                             xy_x *= sinc * this->m_proj_parm.cosph0;
 156                         sinchk:
 157                             if (fabs(lp_lat) >= 1.)
 158                                 lp_lat = lp_lat < 0. ? -half_pi : half_pi;
 159                             else
 160                                 lp_lat = asin(lp_lat);
 161                             break;
 162                         }
 163                         lp_lon = (xy_y == 0. && (this->m_proj_parm.mode == obliq || this->m_proj_parm.mode == equit))
 164                              ? (xy_x == 0. ? 0. : xy_x < 0. ? -half_pi : half_pi)
 165                                            : atan2(xy_x, xy_y);
 166                     }
 167                 }
 168
 169                 static inline std::string get_name()
 170                 {
 171                     return "ortho_spheroid";
 172                 }
 173
 174             };
 175
 176             // Orthographic
 177             template <typename Parameters, typename T>
 178             inline void setup_ortho(Parameters& par, par_ortho<T>& proj_parm)
 179             {
 180                 if (fabs(fabs(par.phi0) - geometry::math::half_pi<T>()) <= epsilon10)
 181                     proj_parm.mode = par.phi0 < 0. ? s_pole : n_pole;
 182                 else if (fabs(par.phi0) > epsilon10) {
 183                     proj_parm.mode = obliq;
 184                     proj_parm.sinph0 = sin(par.phi0);
 185                     proj_parm.cosph0 = cos(par.phi0);
 186                 } else
 187                     proj_parm.mode = equit;
 188                 par.es = 0.;
 189             }
 190
 191     }} // namespace detail::ortho
 192     #endif // doxygen
 193
 194     /*!
 195         \brief Orthographic projection
 196         \ingroup projections
 197         \tparam Geographic latlong point type
 198         \tparam Cartesian xy point type
 199         \tparam Parameters parameter type
 200         \par Projection characteristics
 201          - Azimuthal
 202          - Spheroid
 203         \par Example
 204         \image html ex_ortho.gif
 205     */
 206     template <typename T, typename Parameters>
 207     struct ortho_spheroid : public detail::ortho::base_ortho_spheroid<T, Parameters>
 208     {
 209         template <typename Params>
 210         inline ortho_spheroid(Params const& , Parameters & par)
 211         {
 212             detail::ortho::setup_ortho(par, this->m_proj_parm);
 213         }
 214     };
 215
 216     #ifndef DOXYGEN_NO_DETAIL
 217     namespace detail
 218     {
 219
 220         // Static projection
 221         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_ortho, ortho_spheroid)
 222
 223         // Factory entry(s)
 224         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(ortho_entry, ortho_spheroid)
 225
 226         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(ortho_init)
 227         {
 228             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(ortho, ortho_entry)
 229         }
 230
 231     } // namespace detail
 232     #endif // doxygen
 233
 234 } // namespace projections
 235
 236 }} // namespace boost::geometry
 237
 238 #endif // BOOST_GEOMETRY_PROJECTIONS_ORTHO_HPP
 239