projections/proj/geos.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Copyright (c) 2004   Gerald I. Evenden
0023 // Copyright (c) 2012   Martin Raspaud
0024
0025 // See also (section 4.4.3.2):
0026 //   http://www.eumetsat.int/en/area4/msg/news/us_doc/cgms_03_26.pdf
0027
0028 // Permission is hereby granted, free of charge, to any person obtaining a
0029 // copy of this software and associated documentation files (the "Software"),
0030 // to deal in the Software without restriction, including without limitation
0031 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0032 // and/or sell copies of the Software, and to permit persons to whom the
0033 // Software is furnished to do so, subject to the following conditions:
0034
0035 // The above copyright notice and this permission notice shall be included
0036 // in all copies or substantial portions of the Software.
0037
0038 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0039 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0040 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0041 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0042 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0043 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0044 // DEALINGS IN THE SOFTWARE.
0045
0046 #ifndef BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP
0047 #define BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP
0048
0049 #include <boost/math/special_functions/hypot.hpp>
0050
0051 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0052 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0053 #include <boost/geometry/srs/projections/impl/projects.hpp>
0054 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0055 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0056
0057 namespace boost { namespace geometry
0058 {
0059
0060 namespace projections
0061 {
0062     #ifndef DOXYGEN_NO_DETAIL
0063     namespace detail { namespace geos
0064     {
0065             template <typename T>
0066             struct par_geos
0067             {
0068                 T           h;
0069                 T           radius_p;
0070                 T           radius_p2;
0071                 T           radius_p_inv2;
0072                 T           radius_g;
0073                 T           radius_g_1;
0074                 T           C;
0075                 bool        flip_axis;
0076             };
0077
0078             template <typename T, typename Parameters>
0079             struct base_geos_ellipsoid
0080             {
0081                 par_geos<T> m_proj_parm;
0082
0083                 // FORWARD(e_forward)  ellipsoid
0084                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0085                 inline void fwd(Parameters const& , T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
0086                 {
0087                     T r, Vx, Vy, Vz, tmp;
0088
0089                     /* Calculation of geocentric latitude. */
0090                     lp_lat = atan (this->m_proj_parm.radius_p2 * tan (lp_lat));
0091
0092                     /* Calculation of the three components of the vector from satellite to
0093                     ** position on earth surface (lon,lat).*/
0094                     r = (this->m_proj_parm.radius_p) / boost::math::hypot(this->m_proj_parm.radius_p * cos (lp_lat), sin (lp_lat));
0095                     Vx = r * cos (lp_lon) * cos (lp_lat);
0096                     Vy = r * sin (lp_lon) * cos (lp_lat);
0097                     Vz = r * sin (lp_lat);
0098
0099                     /* Check visibility. */
0100                     if (((this->m_proj_parm.radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * this->m_proj_parm.radius_p_inv2) < 0.) {
0101                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0102                     }
0103
0104                     /* Calculation based on view angles from satellite. */
0105                     tmp = this->m_proj_parm.radius_g - Vx;
0106
0107                     if(this->m_proj_parm.flip_axis) {
0108                         xy_x = this->m_proj_parm.radius_g_1 * atan (Vy / boost::math::hypot (Vz, tmp));
0109                         xy_y = this->m_proj_parm.radius_g_1 * atan (Vz / tmp);
0110                     } else {
0111                         xy_x = this->m_proj_parm.radius_g_1 * atan (Vy / tmp);
0112                         xy_y = this->m_proj_parm.radius_g_1 * atan (Vz / boost::math::hypot (Vy, tmp));
0113                     }
0114                 }
0115
0116                 // INVERSE(e_inverse)  ellipsoid
0117                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0118                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0119                 {
0120                     T Vx, Vy, Vz, a, b, det, k;
0121
0122                     /* Setting three components of vector from satellite to position.*/
0123                     Vx = -1.0;
0124
0125                     if(this->m_proj_parm.flip_axis) {
0126                         Vz = tan (xy_y / this->m_proj_parm.radius_g_1);
0127                         Vy = tan (xy_x / this->m_proj_parm.radius_g_1) * boost::math::hypot(1.0, Vz);
0128                     } else {
0129                         Vy = tan (xy_x / this->m_proj_parm.radius_g_1);
0130                         Vz = tan (xy_y / this->m_proj_parm.radius_g_1) * boost::math::hypot(1.0, Vy);
0131                     }
0132
0133                     /* Calculation of terms in cubic equation and determinant.*/
0134                     a = Vz / this->m_proj_parm.radius_p;
0135                     a   = Vy * Vy + a * a + Vx * Vx;
0136                     b   = 2 * this->m_proj_parm.radius_g * Vx;
0137                     if ((det = (b * b) - 4 * a * this->m_proj_parm.C) < 0.) {
0138                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0139                     }
0140
0141                     /* Calculation of three components of vector from satellite to position.*/
0142                     k  = (-b - sqrt(det)) / (2. * a);
0143                     Vx = this->m_proj_parm.radius_g + k * Vx;
0144                     Vy *= k;
0145                     Vz *= k;
0146
0147                     /* Calculation of longitude and latitude.*/
0148                     lp_lon = atan2 (Vy, Vx);
0149                     lp_lat = atan (Vz * cos (lp_lon) / Vx);
0150                     lp_lat = atan (this->m_proj_parm.radius_p_inv2 * tan (lp_lat));
0151                 }
0152
0153                 static inline std::string get_name()
0154                 {
0155                     return "geos_ellipsoid";
0156                 }
0157
0158             };
0159
0160             template <typename T, typename Parameters>
0161             struct base_geos_spheroid
0162             {
0163                 par_geos<T> m_proj_parm;
0164
0165                 // FORWARD(s_forward)  spheroid
0166                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0167                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0168                 {
0169                     T Vx, Vy, Vz, tmp;
0170
0171                     /* Calculation of the three components of the vector from satellite to
0172                     ** position on earth surface (lon,lat).*/
0173                     tmp = cos(lp_lat);
0174                     Vx = cos (lp_lon) * tmp;
0175                     Vy = sin (lp_lon) * tmp;
0176                     Vz = sin (lp_lat);
0177
0178                     /* Check visibility.*/
0179                     // TODO: in proj4 5.0.0 this check is not present
0180                     if (((this->m_proj_parm.radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.)
0181                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0182
0183                     /* Calculation based on view angles from satellite.*/
0184                     tmp = this->m_proj_parm.radius_g - Vx;
0185
0186                     if(this->m_proj_parm.flip_axis) {
0187                         xy_x = this->m_proj_parm.radius_g_1 * atan(Vy / boost::math::hypot(Vz, tmp));
0188                         xy_y = this->m_proj_parm.radius_g_1 * atan(Vz / tmp);
0189                     } else {
0190                         xy_x = this->m_proj_parm.radius_g_1 * atan(Vy / tmp);
0191                         xy_y = this->m_proj_parm.radius_g_1 * atan(Vz / boost::math::hypot(Vy, tmp));
0192                     }
0193                 }
0194
0195                 // INVERSE(s_inverse)  spheroid
0196                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0197                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0198                 {
0199                     T Vx, Vy, Vz, a, b, det, k;
0200
0201                     /* Setting three components of vector from satellite to position.*/
0202                     Vx = -1.0;
0203                     if(this->m_proj_parm.flip_axis) {
0204                         Vz = tan (xy_y / (this->m_proj_parm.radius_g - 1.0));
0205                         Vy = tan (xy_x / (this->m_proj_parm.radius_g - 1.0)) * sqrt (1.0 + Vz * Vz);
0206                     } else {
0207                         Vy = tan (xy_x / (this->m_proj_parm.radius_g - 1.0));
0208                         Vz = tan (xy_y / (this->m_proj_parm.radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
0209                     }
0210
0211                     /* Calculation of terms in cubic equation and determinant.*/
0212                     a   = Vy * Vy + Vz * Vz + Vx * Vx;
0213                     b   = 2 * this->m_proj_parm.radius_g * Vx;
0214                     if ((det = (b * b) - 4 * a * this->m_proj_parm.C) < 0.) {
0215                         BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
0216                     }
0217
0218                     /* Calculation of three components of vector from satellite to position.*/
0219                     k  = (-b - sqrt(det)) / (2 * a);
0220                     Vx = this->m_proj_parm.radius_g + k * Vx;
0221                     Vy *= k;
0222                     Vz *= k;
0223
0224                     /* Calculation of longitude and latitude.*/
0225                     lp_lon = atan2 (Vy, Vx);
0226                     lp_lat = atan (Vz * cos (lp_lon) / Vx);
0227                 }
0228
0229                 static inline std::string get_name()
0230                 {
0231                     return "geos_spheroid";
0232                 }
0233
0234             };
0235
0236             inline bool geos_flip_axis(srs::detail::proj4_parameters const& params)
0237             {
0238                 std::string sweep_axis = pj_get_param_s(params, "sweep");
0239                 if (sweep_axis.empty())
0240                     return false;
0241                 else {
0242                     if (sweep_axis[1] != '\0' || (sweep_axis[0] != 'x' && sweep_axis[0] != 'y'))
0243                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_sweep_axis) );
0244
0245                     if (sweep_axis[0] == 'x')
0246                         return true;
0247                     else
0248                         return false;
0249                 }
0250             }
0251
0252             template <typename T>
0253             inline bool geos_flip_axis(srs::dpar::parameters<T> const& params)
0254             {
0255                 typename srs::dpar::parameters<T>::const_iterator
0256                     it = pj_param_find(params, srs::dpar::sweep);
0257                 if (it == params.end()) {
0258                     return false;
0259                 } else {
0260                     srs::dpar::value_sweep s = static_cast<srs::dpar::value_sweep>(it->template get_value<int>());
0261                     return s == srs::dpar::sweep_x;
0262                 }
0263             }
0264
0265             // Geostationary Satellite View
0266             template <typename Params, typename Parameters, typename T>
0267             inline void setup_geos(Params const& params, Parameters& par, par_geos<T>& proj_parm)
0268             {
0269                 std::string sweep_axis;
0270
0271                 if ((proj_parm.h = pj_get_param_f<T, srs::spar::h>(params, "h", srs::dpar::h)) <= 0.)
0272                     BOOST_THROW_EXCEPTION( projection_exception(error_h_less_than_zero) );
0273
0274                 if (par.phi0 != 0.0)
0275                     BOOST_THROW_EXCEPTION( projection_exception(error_unknown_prime_meridian) );
0276
0277
0278                 proj_parm.flip_axis = geos_flip_axis(params);
0279
0280                 proj_parm.radius_g_1 = proj_parm.h / par.a;
0281                 proj_parm.radius_g = 1. + proj_parm.radius_g_1;
0282                 proj_parm.C  = proj_parm.radius_g * proj_parm.radius_g - 1.0;
0283                 if (par.es != 0.0) {
0284                     proj_parm.radius_p      = sqrt (par.one_es);
0285                     proj_parm.radius_p2     = par.one_es;
0286                     proj_parm.radius_p_inv2 = par.rone_es;
0287                 } else {
0288                     proj_parm.radius_p = proj_parm.radius_p2 = proj_parm.radius_p_inv2 = 1.0;
0289                 }
0290             }
0291
0292     }} // namespace detail::geos
0293     #endif // doxygen
0294
0295     /*!
0296         \brief Geostationary Satellite View projection
0297         \ingroup projections
0298         \tparam Geographic latlong point type
0299         \tparam Cartesian xy point type
0300         \tparam Parameters parameter type
0301         \par Projection characteristics
0302          - Azimuthal
0303          - Spheroid
0304          - Ellipsoid
0305         \par Projection parameters
0306          - h: Height (real)
0307          - sweep: Sweep axis ('x' or 'y') (string)
0308         \par Example
0309         \image html ex_geos.gif
0310     */
0311     template <typename T, typename Parameters>
0312     struct geos_ellipsoid : public detail::geos::base_geos_ellipsoid<T, Parameters>
0313     {
0314         template <typename Params>
0315         inline geos_ellipsoid(Params const& params, Parameters const& par)
0316         {
0317             detail::geos::setup_geos(params, par, this->m_proj_parm);
0318         }
0319     };
0320
0321     /*!
0322         \brief Geostationary Satellite View projection
0323         \ingroup projections
0324         \tparam Geographic latlong point type
0325         \tparam Cartesian xy point type
0326         \tparam Parameters parameter type
0327         \par Projection characteristics
0328          - Azimuthal
0329          - Spheroid
0330          - Ellipsoid
0331         \par Projection parameters
0332          - h: Height (real)
0333          - sweep: Sweep axis ('x' or 'y') (string)
0334         \par Example
0335         \image html ex_geos.gif
0336     */
0337     template <typename T, typename Parameters>
0338     struct geos_spheroid : public detail::geos::base_geos_spheroid<T, Parameters>
0339     {
0340         template <typename Params>
0341         inline geos_spheroid(Params const& params, Parameters const& par)
0342         {
0343             detail::geos::setup_geos(params, par, this->m_proj_parm);
0344         }
0345     };
0346
0347     #ifndef DOXYGEN_NO_DETAIL
0348     namespace detail
0349     {
0350
0351         // Static projection
0352         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_geos, geos_spheroid, geos_ellipsoid)
0353
0354         // Factory entry(s)
0355         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(geos_entry, geos_spheroid, geos_ellipsoid)
0356
0357         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(geos_init)
0358         {
0359             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(geos, geos_entry);
0360         }
0361
0362     } // namespace detail
0363     #endif // doxygen
0364
0365 } // namespace projections
0366
0367 }} // namespace boost::geometry
0368
0369 #endif // BOOST_GEOMETRY_PROJECTIONS_GEOS_HPP
0370