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) fwdboost::geometry::projections::detail::ortho::base_ortho_spheroid85 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) invboost::geometry::projections::detail::ortho::base_ortho_spheroid122 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 get_nameboost::geometry::projections::detail::ortho::base_ortho_spheroid169 static inline std::string get_name() 170 { 171 return "ortho_spheroid"; 172 } 173 174 }; 175 176 // Orthographic 177 template <typename Parameters, typename T> setup_ortho(Parameters & par,par_ortho<T> & proj_parm)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> ortho_spheroidboost::geometry::projections::ortho_spheroid210 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 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_ortho,ortho_spheroid)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 240