projections/proj/krovak.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, 2022.
0006 // Modifications copyright (c) 2017-2022, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0008 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0009
0010 // Use, modification and distribution is subject to the Boost Software License,
0011 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0012 // http://www.boost.org/LICENSE_1_0.txt)
0013
0014 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0015 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0016 // PROJ4 is maintained by Frank Warmerdam
0017 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0018
0019 // Last updated version of proj: 5.0.0
0020
0021 // Original copyright notice:
0022
0023 // Purpose:  Implementation of the krovak (Krovak) projection.
0024 //           Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3
0025 // Author:   Thomas Flemming, tf@ttqv.com
0026 // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
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_KROVAK_HPP
0047 #define BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP
0048
0049 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0050 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0051 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0052 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0053 #include <boost/geometry/srs/projections/impl/projects.hpp>
0054
0055 namespace boost { namespace geometry
0056 {
0057
0058 namespace projections
0059 {
0060     #ifndef DOXYGEN_NO_DETAIL
0061     namespace detail { namespace krovak
0062     {
0063             static double epsilon = 1e-15;
0064             static double S45 = 0.785398163397448;  /* 45 deg */
0065             static double S90 = 1.570796326794896;  /* 90 deg */
0066             static double UQ  = 1.04216856380474;   /* DU(2, 59, 42, 42.69689) */
0067             static double S0  = 1.37008346281555;   /* Latitude of pseudo standard parallel 78deg 30'00" N */
0068             /* Not sure at all of the appropriate number for max_iter... */
0069             static int max_iter = 100;
0070
0071             template <typename T>
0072             struct par_krovak
0073             {
0074                 T alpha;
0075                 T k;
0076                 T n;
0077                 T rho0;
0078                 T ad;
0079                 int czech;
0080             };
0081
0082             /**
0083                NOTES: According to EPSG the full Krovak projection method should have
0084                       the following parameters.  Within PROJ.4 the azimuth, and pseudo
0085                       standard parallel are hardcoded in the algorithm and can't be
0086                       altered from outside.  The others all have defaults to match the
0087                       common usage with Krovak projection.
0088
0089               lat_0 = latitude of centre of the projection
0090
0091               lon_0 = longitude of centre of the projection
0092
0093               ** = azimuth (true) of the centre line passing through the centre of the projection
0094
0095               ** = latitude of pseudo standard parallel
0096
0097               k  = scale factor on the pseudo standard parallel
0098
0099               x_0 = False Easting of the centre of the projection at the apex of the cone
0100
0101               y_0 = False Northing of the centre of the projection at the apex of the cone
0102
0103              **/
0104
0105             template <typename T, typename Parameters>
0106             struct base_krovak_ellipsoid
0107             {
0108                 par_krovak<T> m_proj_parm;
0109
0110                 // FORWARD(e_forward)  ellipsoid
0111                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0112                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0113                 {
0114                     T gfi, u, deltav, s, d, eps, rho;
0115
0116                     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));
0117
0118                     u = 2. * (atan(this->m_proj_parm.k * math::pow( tan(lp_lat / T(2) + S45), this->m_proj_parm.alpha) / gfi)-S45);
0119                     deltav = -lp_lon * this->m_proj_parm.alpha;
0120
0121                     s = asin(cos(this->m_proj_parm.ad) * sin(u) + sin(this->m_proj_parm.ad) * cos(u) * cos(deltav));
0122                     d = asin(cos(u) * sin(deltav) / cos(s));
0123
0124                     eps = this->m_proj_parm.n * d;
0125                     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);
0126
0127                     xy_y = rho * cos(eps);
0128                     xy_x = rho * sin(eps);
0129
0130                     xy_y *= this->m_proj_parm.czech;
0131                     xy_x *= this->m_proj_parm.czech;
0132                     if (this->m_proj_parm.czech == 1) std::swap(xy_x, xy_y);
0133                 }
0134
0135                 // INVERSE(e_inverse)  ellipsoid
0136                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0137                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0138                 {
0139                     T u, deltav, s, d, eps, rho, fi1;
0140                     int i;
0141
0142                     if (this->m_proj_parm.czech == -1) std::swap(xy_x, xy_y);
0143
0144                     xy_x *= this->m_proj_parm.czech;
0145                     xy_y *= this->m_proj_parm.czech;
0146
0147                     rho = sqrt(xy_x * xy_x + xy_y * xy_y);
0148                     eps = atan2(xy_y, xy_x);
0149
0150                     d = eps / sin(S0);
0151                     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);
0152
0153                     u = asin(cos(this->m_proj_parm.ad) * sin(s) - sin(this->m_proj_parm.ad) * cos(s) * cos(d));
0154                     deltav = asin(cos(s) * sin(d) / cos(u));
0155
0156                     lp_lon = par.lam0 - deltav / this->m_proj_parm.alpha;
0157
0158                     /* ITERATION FOR lp_lat */
0159                     fi1 = u;
0160
0161                     for (i = max_iter; i ; --i) {
0162                         lp_lat = T(2) * ( atan( math::pow( this->m_proj_parm.k, T(-1) / this->m_proj_parm.alpha)  *
0163                                               math::pow( tan(u / T(2) + S45) , T(1) / this->m_proj_parm.alpha)  *
0164                                               math::pow( (T(1) + par.e * sin(fi1)) / (T(1) - par.e * sin(fi1)) , par.e / T(2))
0165                                             )  - S45);
0166
0167                         if (fabs(fi1 - lp_lat) < epsilon)
0168                             break;
0169                         fi1 = lp_lat;
0170                     }
0171                     if( i == 0 )
0172                         BOOST_THROW_EXCEPTION( projection_exception(error_non_convergent) );
0173
0174                    lp_lon -= par.lam0;
0175                 }
0176
0177                 static inline std::string get_name()
0178                 {
0179                     return "krovak_ellipsoid";
0180                 }
0181
0182             };
0183
0184             // Krovak
0185             template <typename Params, typename Parameters, typename T>
0186             inline void setup_krovak(Params const& params, Parameters& par, par_krovak<T>& proj_parm)
0187             {
0188                 T u0, n0, g;
0189
0190                 /* we want Bessel as fixed ellipsoid */
0191                 par.a = 6377397.155;
0192                 par.es = 0.006674372230614;
0193                 par.e = sqrt(par.es);
0194
0195                 /* if latitude of projection center is not set, use 49d30'N */
0196                 if (!pj_param_exists<srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0))
0197                     par.phi0 = 0.863937979737193;
0198
0199                 /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
0200                 /* that will correspond to using longitudes relative to greenwich    */
0201                 /* as input and output, instead of lat/long relative to Ferro */
0202                 if (!pj_param_exists<srs::spar::lon_0>(params, "lon_0", srs::dpar::lon_0))
0203                     par.lam0 = 0.7417649320975901 - 0.308341501185665;
0204
0205                 /* if scale not set default to 0.9999 */
0206                 if (!pj_param_exists<srs::spar::k>(params, "k", srs::dpar::k))
0207                     par.k0 = 0.9999;
0208
0209                 proj_parm.czech = 1;
0210                 if( !pj_param_exists<srs::spar::czech>(params, "czech", srs::dpar::czech) )
0211                     proj_parm.czech = -1;
0212
0213                 /* Set up shared parameters between forward and inverse */
0214                 proj_parm.alpha = sqrt(T(1) + (par.es * math::pow(cos(par.phi0), 4)) / (T(1) - par.es));
0215                 u0 = asin(sin(par.phi0) / proj_parm.alpha);
0216                 g = math::pow( (T(1) + par.e * sin(par.phi0)) / (T(1) - par.e * sin(par.phi0)) , proj_parm.alpha * par.e / T(2) );
0217                 proj_parm.k = tan( u0 / 2. + S45) / math::pow(tan(par.phi0 / T(2) + S45) , proj_parm.alpha) * g;
0218                 n0 = sqrt(T(1) - par.es) / (T(1) - par.es * math::pow(sin(par.phi0), 2));
0219                 proj_parm.n = sin(S0);
0220                 proj_parm.rho0 = par.k0 * n0 / tan(S0);
0221                 proj_parm.ad = S90 - UQ;
0222             }
0223
0224     }} // namespace detail::krovak
0225     #endif // doxygen
0226
0227     /*!
0228         \brief Krovak projection
0229         \ingroup projections
0230         \tparam Geographic latlong point type
0231         \tparam Cartesian xy point type
0232         \tparam Parameters parameter type
0233         \par Projection characteristics
0234          - Pseudocylindrical
0235          - Ellipsoid
0236         \par Projection parameters
0237          - lat_ts: Latitude of true scale (degrees)
0238          - lat_0: Latitude of origin
0239          - lon_0: Central meridian
0240          - k: Scale factor on the pseudo standard parallel
0241         \par Example
0242         \image html ex_krovak.gif
0243     */
0244     template <typename T, typename Parameters>
0245     struct krovak_ellipsoid : public detail::krovak::base_krovak_ellipsoid<T, Parameters>
0246     {
0247         template <typename Params>
0248         inline krovak_ellipsoid(Params const& params, Parameters & par)
0249         {
0250             detail::krovak::setup_krovak(params, par, this->m_proj_parm);
0251         }
0252     };
0253
0254     #ifndef DOXYGEN_NO_DETAIL
0255     namespace detail
0256     {
0257
0258         // Static projection
0259         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_krovak, krovak_ellipsoid)
0260
0261         // Factory entry(s)
0262         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(krovak_entry, krovak_ellipsoid)
0263
0264         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(krovak_init)
0265         {
0266             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(krovak, krovak_entry)
0267         }
0268
0269     } // namespace detail
0270     #endif // doxygen
0271
0272 } // namespace projections
0273
0274 }} // namespace boost::geometry
0275
0276 #endif // BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP
0277