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 // Copyright (c) 2004 Gerald I. Evenden 23 // Copyright (c) 2012 Martin Raspaud 24 25 // See also (section 4.4.3.2): 26 // http://www.eumetsat.int/en/area4/msg/news/us_doc/cgms_03_26.pdf 27 28 // Permission is hereby granted, free of charge, to any person obtaining a 29 // copy of this software and associated documentation files (the "Software"), 30 // to deal in the Software without restriction, including without limitation 31 // the rights to use, copy, modify, merge, publish, distribute, sublicense, 32 // and/or sell copies of the Software, and to permit persons to whom the 33 // Software is furnished to do so, subject to the following conditions: 34 35 // The above copyright notice and this permission notice shall be included 36 // in all copies or substantial portions of the Software. 37 38 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 39 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 40 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 41 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 42 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 43 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 44 // DEALINGS IN THE SOFTWARE. 45 46 #ifndef BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP 47 #define BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP 48 49 #include <boost/math/special_functions/hypot.hpp> 50 51 #include <boost/geometry/srs/projections/impl/base_static.hpp> 52 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> 53 #include <boost/geometry/srs/projections/impl/projects.hpp> 54 #include <boost/geometry/srs/projections/impl/factory_entry.hpp> 55 #include <boost/geometry/srs/projections/impl/pj_param.hpp> 56 57 namespace boost { namespace geometry 58 { 59 60 namespace projections 61 { 62 #ifndef DOXYGEN_NO_DETAIL 63 namespace detail { namespace geos 64 { 65 template <typename T> 66 struct par_geos 67 { 68 T h; 69 T radius_p; 70 T radius_p2; 71 T radius_p_inv2; 72 T radius_g; 73 T radius_g_1; 74 T C; 75 bool flip_axis; 76 }; 77 78 template <typename T, typename Parameters> 79 struct base_geos_ellipsoid 80 { 81 par_geos<T> m_proj_parm; 82 83 // FORWARD(e_forward) ellipsoid 84 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::geos::base_geos_ellipsoid85 inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const 86 { 87 T r, Vx, Vy, Vz, tmp; 88 89 /* Calculation of geocentric latitude. */ 90 lp_lat = atan (this->m_proj_parm.radius_p2 * tan (lp_lat)); 91 92 /* Calculation of the three components of the vector from satellite to 93 ** position on earth surface (lon,lat).*/ 94 r = (this->m_proj_parm.radius_p) / boost::math::hypot(this->m_proj_parm.radius_p * cos (lp_lat), sin (lp_lat)); 95 Vx = r * cos (lp_lon) * cos (lp_lat); 96 Vy = r * sin (lp_lon) * cos (lp_lat); 97 Vz = r * sin (lp_lat); 98 99 /* Check visibility. */ 100 if (((this->m_proj_parm.radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * this->m_proj_parm.radius_p_inv2) < 0.) { 101 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 102 } 103 104 /* Calculation based on view angles from satellite. */ 105 tmp = this->m_proj_parm.radius_g - Vx; 106 107 if(this->m_proj_parm.flip_axis) { 108 xy_x = this->m_proj_parm.radius_g_1 * atan (Vy / boost::math::hypot (Vz, tmp)); 109 xy_y = this->m_proj_parm.radius_g_1 * atan (Vz / tmp); 110 } else { 111 xy_x = this->m_proj_parm.radius_g_1 * atan (Vy / tmp); 112 xy_y = this->m_proj_parm.radius_g_1 * atan (Vz / boost::math::hypot (Vy, tmp)); 113 } 114 } 115 116 // INVERSE(e_inverse) ellipsoid 117 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::geos::base_geos_ellipsoid118 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 119 { 120 T Vx, Vy, Vz, a, b, det, k; 121 122 /* Setting three components of vector from satellite to position.*/ 123 Vx = -1.0; 124 125 if(this->m_proj_parm.flip_axis) { 126 Vz = tan (xy_y / this->m_proj_parm.radius_g_1); 127 Vy = tan (xy_x / this->m_proj_parm.radius_g_1) * boost::math::hypot(1.0, Vz); 128 } else { 129 Vy = tan (xy_x / this->m_proj_parm.radius_g_1); 130 Vz = tan (xy_y / this->m_proj_parm.radius_g_1) * boost::math::hypot(1.0, Vy); 131 } 132 133 /* Calculation of terms in cubic equation and determinant.*/ 134 a = Vz / this->m_proj_parm.radius_p; 135 a = Vy * Vy + a * a + Vx * Vx; 136 b = 2 * this->m_proj_parm.radius_g * Vx; 137 if ((det = (b * b) - 4 * a * this->m_proj_parm.C) < 0.) { 138 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 139 } 140 141 /* Calculation of three components of vector from satellite to position.*/ 142 k = (-b - sqrt(det)) / (2. * a); 143 Vx = this->m_proj_parm.radius_g + k * Vx; 144 Vy *= k; 145 Vz *= k; 146 147 /* Calculation of longitude and latitude.*/ 148 lp_lon = atan2 (Vy, Vx); 149 lp_lat = atan (Vz * cos (lp_lon) / Vx); 150 lp_lat = atan (this->m_proj_parm.radius_p_inv2 * tan (lp_lat)); 151 } 152 get_nameboost::geometry::projections::detail::geos::base_geos_ellipsoid153 static inline std::string get_name() 154 { 155 return "geos_ellipsoid"; 156 } 157 158 }; 159 160 template <typename T, typename Parameters> 161 struct base_geos_spheroid 162 { 163 par_geos<T> m_proj_parm; 164 165 // FORWARD(s_forward) spheroid 166 // Project coordinates from geographic (lon, lat) to cartesian (x, y) fwdboost::geometry::projections::detail::geos::base_geos_spheroid167 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const 168 { 169 T Vx, Vy, Vz, tmp; 170 171 /* Calculation of the three components of the vector from satellite to 172 ** position on earth surface (lon,lat).*/ 173 tmp = cos(lp_lat); 174 Vx = cos (lp_lon) * tmp; 175 Vy = sin (lp_lon) * tmp; 176 Vz = sin (lp_lat); 177 178 /* Check visibility.*/ 179 // TODO: in proj4 5.0.0 this check is not present 180 if (((this->m_proj_parm.radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) 181 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 182 183 /* Calculation based on view angles from satellite.*/ 184 tmp = this->m_proj_parm.radius_g - Vx; 185 186 if(this->m_proj_parm.flip_axis) { 187 xy_x = this->m_proj_parm.radius_g_1 * atan(Vy / boost::math::hypot(Vz, tmp)); 188 xy_y = this->m_proj_parm.radius_g_1 * atan(Vz / tmp); 189 } else { 190 xy_x = this->m_proj_parm.radius_g_1 * atan(Vy / tmp); 191 xy_y = this->m_proj_parm.radius_g_1 * atan(Vz / boost::math::hypot(Vy, tmp)); 192 } 193 } 194 195 // INVERSE(s_inverse) spheroid 196 // Project coordinates from cartesian (x, y) to geographic (lon, lat) invboost::geometry::projections::detail::geos::base_geos_spheroid197 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const 198 { 199 T Vx, Vy, Vz, a, b, det, k; 200 201 /* Setting three components of vector from satellite to position.*/ 202 Vx = -1.0; 203 if(this->m_proj_parm.flip_axis) { 204 Vz = tan (xy_y / (this->m_proj_parm.radius_g - 1.0)); 205 Vy = tan (xy_x / (this->m_proj_parm.radius_g - 1.0)) * sqrt (1.0 + Vz * Vz); 206 } else { 207 Vy = tan (xy_x / (this->m_proj_parm.radius_g - 1.0)); 208 Vz = tan (xy_y / (this->m_proj_parm.radius_g - 1.0)) * sqrt (1.0 + Vy * Vy); 209 } 210 211 /* Calculation of terms in cubic equation and determinant.*/ 212 a = Vy * Vy + Vz * Vz + Vx * Vx; 213 b = 2 * this->m_proj_parm.radius_g * Vx; 214 if ((det = (b * b) - 4 * a * this->m_proj_parm.C) < 0.) { 215 BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) ); 216 } 217 218 /* Calculation of three components of vector from satellite to position.*/ 219 k = (-b - sqrt(det)) / (2 * a); 220 Vx = this->m_proj_parm.radius_g + k * Vx; 221 Vy *= k; 222 Vz *= k; 223 224 /* Calculation of longitude and latitude.*/ 225 lp_lon = atan2 (Vy, Vx); 226 lp_lat = atan (Vz * cos (lp_lon) / Vx); 227 } 228 get_nameboost::geometry::projections::detail::geos::base_geos_spheroid229 static inline std::string get_name() 230 { 231 return "geos_spheroid"; 232 } 233 234 }; 235 geos_flip_axis(srs::detail::proj4_parameters const & params)236 inline bool geos_flip_axis(srs::detail::proj4_parameters const& params) 237 { 238 std::string sweep_axis = pj_get_param_s(params, "sweep"); 239 if (sweep_axis.empty()) 240 return false; 241 else { 242 if (sweep_axis[1] != '\0' || (sweep_axis[0] != 'x' && sweep_axis[0] != 'y')) 243 BOOST_THROW_EXCEPTION( projection_exception(error_invalid_sweep_axis) ); 244 245 if (sweep_axis[0] == 'x') 246 return true; 247 else 248 return false; 249 } 250 } 251 252 template <typename T> geos_flip_axis(srs::dpar::parameters<T> const & params)253 inline bool geos_flip_axis(srs::dpar::parameters<T> const& params) 254 { 255 typename srs::dpar::parameters<T>::const_iterator 256 it = pj_param_find(params, srs::dpar::sweep); 257 if (it == params.end()) { 258 return false; 259 } else { 260 srs::dpar::value_sweep s = static_cast<srs::dpar::value_sweep>(it->template get_value<int>()); 261 return s == srs::dpar::sweep_x; 262 } 263 } 264 265 // Geostationary Satellite View 266 template <typename Params, typename Parameters, typename T> setup_geos(Params const & params,Parameters & par,par_geos<T> & proj_parm)267 inline void setup_geos(Params const& params, Parameters& par, par_geos<T>& proj_parm) 268 { 269 std::string sweep_axis; 270 271 if ((proj_parm.h = pj_get_param_f<T, srs::spar::h>(params, "h", srs::dpar::h)) <= 0.) 272 BOOST_THROW_EXCEPTION( projection_exception(error_h_less_than_zero) ); 273 274 if (par.phi0 != 0.0) 275 BOOST_THROW_EXCEPTION( projection_exception(error_unknown_prime_meridian) ); 276 277 278 proj_parm.flip_axis = geos_flip_axis(params); 279 280 proj_parm.radius_g_1 = proj_parm.h / par.a; 281 proj_parm.radius_g = 1. + proj_parm.radius_g_1; 282 proj_parm.C = proj_parm.radius_g * proj_parm.radius_g - 1.0; 283 if (par.es != 0.0) { 284 proj_parm.radius_p = sqrt (par.one_es); 285 proj_parm.radius_p2 = par.one_es; 286 proj_parm.radius_p_inv2 = par.rone_es; 287 } else { 288 proj_parm.radius_p = proj_parm.radius_p2 = proj_parm.radius_p_inv2 = 1.0; 289 } 290 } 291 292 }} // namespace detail::geos 293 #endif // doxygen 294 295 /*! 296 \brief Geostationary Satellite View projection 297 \ingroup projections 298 \tparam Geographic latlong point type 299 \tparam Cartesian xy point type 300 \tparam Parameters parameter type 301 \par Projection characteristics 302 - Azimuthal 303 - Spheroid 304 - Ellipsoid 305 \par Projection parameters 306 - h: Height (real) 307 - sweep: Sweep axis ('x' or 'y') (string) 308 \par Example 309 \image html ex_geos.gif 310 */ 311 template <typename T, typename Parameters> 312 struct geos_ellipsoid : public detail::geos::base_geos_ellipsoid<T, Parameters> 313 { 314 template <typename Params> geos_ellipsoidboost::geometry::projections::geos_ellipsoid315 inline geos_ellipsoid(Params const& params, Parameters const& par) 316 { 317 detail::geos::setup_geos(params, par, this->m_proj_parm); 318 } 319 }; 320 321 /*! 322 \brief Geostationary Satellite View projection 323 \ingroup projections 324 \tparam Geographic latlong point type 325 \tparam Cartesian xy point type 326 \tparam Parameters parameter type 327 \par Projection characteristics 328 - Azimuthal 329 - Spheroid 330 - Ellipsoid 331 \par Projection parameters 332 - h: Height (real) 333 - sweep: Sweep axis ('x' or 'y') (string) 334 \par Example 335 \image html ex_geos.gif 336 */ 337 template <typename T, typename Parameters> 338 struct geos_spheroid : public detail::geos::base_geos_spheroid<T, Parameters> 339 { 340 template <typename Params> geos_spheroidboost::geometry::projections::geos_spheroid341 inline geos_spheroid(Params const& params, Parameters const& par) 342 { 343 detail::geos::setup_geos(params, par, this->m_proj_parm); 344 } 345 }; 346 347 #ifndef DOXYGEN_NO_DETAIL 348 namespace detail 349 { 350 351 // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_geos,geos_spheroid,geos_ellipsoid)352 BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_geos, geos_spheroid, geos_ellipsoid) 353 354 // Factory entry(s) 355 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(geos_entry, geos_spheroid, geos_ellipsoid) 356 357 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(geos_init) 358 { 359 BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(geos, geos_entry); 360 } 361 362 } // namespace detail 363 #endif // doxygen 364 365 } // namespace projections 366 367 }} // namespace boost::geometry 368 369 #endif // BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP 370 371