boost/geometry/srs/projections/proj/krovak.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 // Purpose:  Implementation of the krovak (Krovak) projection.
  23 //           Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3
  24 // Author:   Thomas Flemming, tf@ttqv.com
  25 // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
  26
  27 // Permission is hereby granted, free of charge, to any person obtaining a
  28 // copy of this software and associated documentation files (the "Software"),
  29 // to deal in the Software without restriction, including without limitation
  30 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  31 // and/or sell copies of the Software, and to permit persons to whom the
  32 // Software is furnished to do so, subject to the following conditions:
  33
  34 // The above copyright notice and this permission notice shall be included
  35 // in all copies or substantial portions of the Software.
  36
  37 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  38 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  39 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  40 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  41 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  42 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  43 // DEALINGS IN THE SOFTWARE.
  44
  45 #ifndef BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP
  46 #define BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP
  47
  48 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  49 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  50 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  51 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  52 #include <boost/geometry/srs/projections/impl/projects.hpp>
  53
  54 namespace boost { namespace geometry
  55 {
  56
  57 namespace projections
  58 {
  59     #ifndef DOXYGEN_NO_DETAIL
  60     namespace detail { namespace krovak
  61     {
  62             static double epsilon = 1e-15;
  63             static double S45 = 0.785398163397448;  /* 45 deg */
  64             static double S90 = 1.570796326794896;  /* 90 deg */
  65             static double UQ  = 1.04216856380474;   /* DU(2, 59, 42, 42.69689) */
  66             static double S0  = 1.37008346281555;   /* Latitude of pseudo standard parallel 78deg 30'00" N */
  67             /* Not sure at all of the appropriate number for max_iter... */
  68             static int max_iter = 100;
  69
  70             template <typename T>
  71             struct par_krovak
  72             {
  73                 T alpha;
  74                 T k;
  75                 T n;
  76                 T rho0;
  77                 T ad;
  78                 int czech;
  79             };
  80
  81             /**
  82                NOTES: According to EPSG the full Krovak projection method should have
  83                       the following parameters.  Within PROJ.4 the azimuth, and pseudo
  84                       standard parallel are hardcoded in the algorithm and can't be
  85                       altered from outside.  The others all have defaults to match the
  86                       common usage with Krovak projection.
  87
  88               lat_0 = latitude of centre of the projection
  89
  90               lon_0 = longitude of centre of the projection
  91
  92               ** = azimuth (true) of the centre line passing through the centre of the projection
  93
  94               ** = latitude of pseudo standard parallel
  95
  96               k  = scale factor on the pseudo standard parallel
  97
  98               x_0 = False Easting of the centre of the projection at the apex of the cone
  99
 100               y_0 = False Northing of the centre of the projection at the apex of the cone
 101
 102              **/
 103
 104             template <typename T, typename Parameters>
 105             struct base_krovak_ellipsoid
 106             {
 107                 par_krovak<T> m_proj_parm;
 108
 109                 // FORWARD(e_forward)  ellipsoid
 110                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 111                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 112                 {
 113                     T gfi, u, deltav, s, d, eps, rho;
 114
 115                     gfi = math::pow( (T(1) + par.e * sin(lp_lat)) / (T(1) - par.e * sin(lp_lat)), this->m_proj_parm.alpha * par.e / T(2));
 116
 117                     u = 2. * (atan(this->m_proj_parm.k * math::pow( tan(lp_lat / T(2) + S45), this->m_proj_parm.alpha) / gfi)-S45);
 118                     deltav = -lp_lon * this->m_proj_parm.alpha;
 119
 120                     s = asin(cos(this->m_proj_parm.ad) * sin(u) + sin(this->m_proj_parm.ad) * cos(u) * cos(deltav));
 121                     d = asin(cos(u) * sin(deltav) / cos(s));
 122
 123                     eps = this->m_proj_parm.n * d;
 124                     rho = this->m_proj_parm.rho0 * math::pow(tan(S0 / T(2) + S45) , this->m_proj_parm.n) / math::pow(tan(s / T(2) + S45) , this->m_proj_parm.n);
 125
 126                     xy_y = rho * cos(eps);
 127                     xy_x = rho * sin(eps);
 128
 129                     xy_y *= this->m_proj_parm.czech;
 130                     xy_x *= this->m_proj_parm.czech;
 131                 }
 132
 133                 // INVERSE(e_inverse)  ellipsoid
 134                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 135                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 136                 {
 137                     T u, deltav, s, d, eps, rho, fi1, xy0;
 138                     int i;
 139
 140                     // TODO: replace with std::swap()
 141                     xy0 = xy_x;
 142                     xy_x = xy_y;
 143                     xy_y = xy0;
 144
 145                     xy_x *= this->m_proj_parm.czech;
 146                     xy_y *= this->m_proj_parm.czech;
 147
 148                     rho = sqrt(xy_x * xy_x + xy_y * xy_y);
 149                     eps = atan2(xy_y, xy_x);
 150
 151                     d = eps / sin(S0);
 152                     s = T(2) * (atan(math::pow(this->m_proj_parm.rho0 / rho, T(1) / this->m_proj_parm.n) * tan(S0 / T(2) + S45)) - S45);
 153
 154                     u = asin(cos(this->m_proj_parm.ad) * sin(s) - sin(this->m_proj_parm.ad) * cos(s) * cos(d));
 155                     deltav = asin(cos(s) * sin(d) / cos(u));
 156
 157                     lp_lon = par.lam0 - deltav / this->m_proj_parm.alpha;
 158
 159                     /* ITERATION FOR lp_lat */
 160                     fi1 = u;
 161
 162                     for (i = max_iter; i ; --i) {
 163                         lp_lat = T(2) * ( atan( math::pow( this->m_proj_parm.k, T(-1) / this->m_proj_parm.alpha)  *
 164                                               math::pow( tan(u / T(2) + S45) , T(1) / this->m_proj_parm.alpha)  *
 165                                               math::pow( (T(1) + par.e * sin(fi1)) / (T(1) - par.e * sin(fi1)) , par.e / T(2))
 166                                             )  - S45);
 167
 168                         if (fabs(fi1 - lp_lat) < epsilon)
 169                             break;
 170                         fi1 = lp_lat;
 171                     }
 172                     if( i == 0 )
 173                         BOOST_THROW_EXCEPTION( projection_exception(error_non_convergent) );
 174
 175                    lp_lon -= par.lam0;
 176                 }
 177
 178                 static inline std::string get_name()
 179                 {
 180                     return "krovak_ellipsoid";
 181                 }
 182
 183             };
 184
 185             // Krovak
 186             template <typename Params, typename Parameters, typename T>
 187             inline void setup_krovak(Params const& params, Parameters& par, par_krovak<T>& proj_parm)
 188             {
 189                 T u0, n0, g;
 190
 191                 /* we want Bessel as fixed ellipsoid */
 192                 par.a = 6377397.155;
 193                 par.es = 0.006674372230614;
 194                 par.e = sqrt(par.es);
 195
 196                 /* if latitude of projection center is not set, use 49d30'N */
 197                 if (!pj_param_exists<srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0))
 198                     par.phi0 = 0.863937979737193;
 199
 200                 /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
 201                 /* that will correspond to using longitudes relative to greenwich    */
 202                 /* as input and output, instead of lat/long relative to Ferro */
 203                 if (!pj_param_exists<srs::spar::lon_0>(params, "lon_0", srs::dpar::lon_0))
 204                     par.lam0 = 0.7417649320975901 - 0.308341501185665;
 205
 206                 /* if scale not set default to 0.9999 */
 207                 if (!pj_param_exists<srs::spar::k>(params, "k", srs::dpar::k))
 208                     par.k0 = 0.9999;
 209
 210                 proj_parm.czech = 1;
 211                 if( !pj_param_exists<srs::spar::czech>(params, "czech", srs::dpar::czech) )
 212                     proj_parm.czech = -1;
 213
 214                 /* Set up shared parameters between forward and inverse */
 215                 proj_parm.alpha = sqrt(T(1) + (par.es * math::pow(cos(par.phi0), 4)) / (T(1) - par.es));
 216                 u0 = asin(sin(par.phi0) / proj_parm.alpha);
 217                 g = math::pow( (T(1) + par.e * sin(par.phi0)) / (T(1) - par.e * sin(par.phi0)) , proj_parm.alpha * par.e / T(2) );
 218                 proj_parm.k = tan( u0 / 2. + S45) / math::pow(tan(par.phi0 / T(2) + S45) , proj_parm.alpha) * g;
 219                 n0 = sqrt(T(1) - par.es) / (T(1) - par.es * math::pow(sin(par.phi0), 2));
 220                 proj_parm.n = sin(S0);
 221                 proj_parm.rho0 = par.k0 * n0 / tan(S0);
 222                 proj_parm.ad = S90 - UQ;
 223             }
 224
 225     }} // namespace detail::krovak
 226     #endif // doxygen
 227
 228     /*!
 229         \brief Krovak projection
 230         \ingroup projections
 231         \tparam Geographic latlong point type
 232         \tparam Cartesian xy point type
 233         \tparam Parameters parameter type
 234         \par Projection characteristics
 235          - Pseudocylindrical
 236          - Ellipsoid
 237         \par Projection parameters
 238          - lat_ts: Latitude of true scale (degrees)
 239          - lat_0: Latitude of origin
 240          - lon_0: Central meridian
 241          - k: Scale factor on the pseudo standard parallel
 242         \par Example
 243         \image html ex_krovak.gif
 244     */
 245     template <typename T, typename Parameters>
 246     struct krovak_ellipsoid : public detail::krovak::base_krovak_ellipsoid<T, Parameters>
 247     {
 248         template <typename Params>
 249         inline krovak_ellipsoid(Params const& params, Parameters & par)
 250         {
 251             detail::krovak::setup_krovak(params, par, this->m_proj_parm);
 252         }
 253     };
 254
 255     #ifndef DOXYGEN_NO_DETAIL
 256     namespace detail
 257     {
 258
 259         // Static projection
 260         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_krovak, krovak_ellipsoid)
 261
 262         // Factory entry(s)
 263         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(krovak_entry, krovak_ellipsoid)
 264
 265         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(krovak_init)
 266         {
 267             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(krovak, krovak_entry)
 268         }
 269
 270     } // namespace detail
 271     #endif // doxygen
 272
 273 } // namespace projections
 274
 275 }} // namespace boost::geometry
 276
 277 #endif // BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP
 278