boost/geometry/srs/projections/proj/geos.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 // 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)
  85                 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)
 118                 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
 153                 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)
 167                 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)
 197                 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
 229                 static inline std::string get_name()
 230                 {
 231                     return "geos_spheroid";
 232                 }
 233
 234             };
 235
 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>
 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>
 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>
 315         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>
 341         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
 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