diff options
Diffstat (limited to 'boost/geometry/srs/projections/proj/krovak.hpp')
-rw-r--r-- | boost/geometry/srs/projections/proj/krovak.hpp | 286 |
1 files changed, 111 insertions, 175 deletions
diff --git a/boost/geometry/srs/projections/proj/krovak.hpp b/boost/geometry/srs/projections/proj/krovak.hpp index 09c24772ed..a008f181db 100644 --- a/boost/geometry/srs/projections/proj/krovak.hpp +++ b/boost/geometry/srs/projections/proj/krovak.hpp @@ -1,13 +1,9 @@ -#ifndef BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP -#define BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP - -// Boost.Geometry - extensions-gis-projections (based on PROJ4) -// This file is automatically generated. DO NOT EDIT. +// Boost.Geometry - gis-projections (based on PROJ4) // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2017. -// Modifications copyright (c) 2017, Oracle and/or its affiliates. +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle. // Use, modification and distribution is subject to the Boost Software License, @@ -19,12 +15,12 @@ // PROJ4 is maintained by Frank Warmerdam // PROJ4 is converted to Boost.Geometry by Barend Gehrels -// Last updated version of proj: 4.9.1 +// Last updated version of proj: 5.0.0 // Original copyright notice: // Purpose: Implementation of the krovak (Krovak) projection. -// Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3 +// Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3 // Author: Thomas Flemming, tf@ttqv.com // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com @@ -46,6 +42,9 @@ // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER // DEALINGS IN THE SOFTWARE. +#ifndef BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP +#define BOOST_GEOMETRY_PROJECTIONS_KROVAK_HPP + #include <boost/geometry/srs/projections/impl/base_static.hpp> #include <boost/geometry/srs/projections/impl/base_dynamic.hpp> #include <boost/geometry/srs/projections/impl/projects.hpp> @@ -56,7 +55,7 @@ namespace boost { namespace geometry namespace srs { namespace par4 { - struct krovak {}; + struct krovak {}; // Krovak }} //namespace srs::par4 @@ -65,10 +64,23 @@ namespace projections #ifndef DOXYGEN_NO_DETAIL namespace detail { namespace krovak { + static double epsilon = 1e-15; + static double S45 = 0.785398163397448; /* 45 deg */ + static double S90 = 1.570796326794896; /* 90 deg */ + static double UQ = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ + static double S0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78deg 30'00" N */ + /* Not sure at all of the appropriate number for max_iter... */ + static int max_iter = 100; + template <typename T> struct par_krovak { - T C_x; + T alpha; + T k; + T n; + T rho0; + T ad; + int czech; }; /** @@ -95,164 +107,81 @@ namespace projections **/ // template class, using CRTP to implement forward/inverse - template <typename CalculationType, typename Parameters> - struct base_krovak_ellipsoid : public base_t_fi<base_krovak_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters> + template <typename T, typename Parameters> + struct base_krovak_ellipsoid + : public base_t_fi<base_krovak_ellipsoid<T, Parameters>, T, Parameters> { - - typedef CalculationType geographic_type; - typedef CalculationType cartesian_type; - - par_krovak<CalculationType> m_proj_parm; + par_krovak<T> m_proj_parm; inline base_krovak_ellipsoid(const Parameters& par) - : base_t_fi<base_krovak_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters>(*this, par) {} + : base_t_fi<base_krovak_ellipsoid<T, Parameters>, T, Parameters>(*this, par) + {} // FORWARD(e_forward) ellipsoid // Project coordinates from geographic (lon, lat) to cartesian (x, y) - inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const + inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const { - /* calculate xy from lat/lon */ - - /* Constants, identical to inverse transform function */ - CalculationType s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n; - CalculationType gfi, u, fi0, deltav, s, d, eps, ro; - - - s45 = 0.785398163397448; /* 45 DEG */ - s90 = 2 * s45; - fi0 = this->m_par.phi0; /* Latitude of projection centre 49 DEG 30' */ - - /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must - be set to 1 here. - Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128, - e2=0.006674372230614; - */ - a = 1; /* 6377397.155; */ - /* e2 = this->m_par.es;*/ /* 0.006674372230614; */ - e2 = 0.006674372230614; - e = sqrt(e2); - - alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2)); + T gfi, u, deltav, s, d, eps, rho; - uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ - u0 = asin(sin(fi0) / alfa); - g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. ); + gfi = math::pow( (T(1) + this->m_par.e * sin(lp_lat)) / (T(1) - this->m_par.e * sin(lp_lat)), this->m_proj_parm.alpha * this->m_par.e / T(2)); - k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g; + u = 2. * (atan(this->m_proj_parm.k * math::pow( tan(lp_lat / T(2) + S45), this->m_proj_parm.alpha) / gfi)-S45); + deltav = -lp_lon * this->m_proj_parm.alpha; - k1 = this->m_par.k0; - n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2)); - s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78 DEG 30'00" N */ - n = sin(s0); - ro0 = k1 * n0 / tan(s0); - ad = s90 - uq; - - /* Transformation */ - - gfi =pow ( ((1. + e * sin(lp_lat)) / - (1. - e * sin(lp_lat))) , (alfa * e / 2.)); + s = asin(cos(this->m_proj_parm.ad) * sin(u) + sin(this->m_proj_parm.ad) * cos(u) * cos(deltav)); + d = asin(cos(u) * sin(deltav) / cos(s)); - u= 2. * (atan(k * pow( tan(lp_lat / 2. + s45), alfa) / gfi)-s45); + eps = this->m_proj_parm.n * d; + 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); - deltav = - lp_lon * alfa; + xy_y = rho * cos(eps); + xy_x = rho * sin(eps); - s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav)); - d = asin(cos(u) * sin(deltav) / cos(s)); - eps = n * d; - ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n) ; - - /* x and y are reverted! */ - xy_y = ro * cos(eps) / a; - xy_x = ro * sin(eps) / a; - - if( !pj_param(this->m_par.params, "tczech").i ) - { - xy_y *= -1.0; - xy_x *= -1.0; - } + xy_y *= this->m_proj_parm.czech; + xy_x *= this->m_proj_parm.czech; } // INVERSE(e_inverse) ellipsoid // Project coordinates from cartesian (x, y) to geographic (lon, lat) - inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const + inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const { - /* calculate lat/lon from xy */ - - /* Constants, identisch wie in der Umkehrfunktion */ - CalculationType s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n; - CalculationType u, deltav, s, d, eps, ro, fi1, xy0; - int ok; - - s45 = 0.785398163397448; /* 45 DEG */ - s90 = 2 * s45; - fi0 = this->m_par.phi0; /* Latitude of projection centre 49 DEG 30' */ - - - /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must - be set to 1 here. - Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128, - e2=0.006674372230614; - */ - a = 1; /* 6377397.155; */ - /* e2 = this->m_par.es; */ /* 0.006674372230614; */ - e2 = 0.006674372230614; - e = sqrt(e2); - - alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2)); - uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ - u0 = asin(sin(fi0) / alfa); - g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. ); - - k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g; - - k1 = this->m_par.k0; - n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2)); - s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78 DEG 30'00" N */ - n = sin(s0); - ro0 = k1 * n0 / tan(s0); - ad = s90 - uq; - - - /* Transformation */ - /* revert y, x*/ - xy0=xy_x; - xy_x=xy_y; - xy_y=xy0; - - if( !pj_param(this->m_par.params, "tczech").i ) - { - xy_x *= -1.0; - xy_y *= -1.0; - } - - ro = sqrt(xy_x * xy_x + xy_y * xy_y); + T u, deltav, s, d, eps, rho, fi1, xy0; + int i; + + // TODO: replace with std::swap() + xy0 = xy_x; + xy_x = xy_y; + xy_y = xy0; + + xy_x *= this->m_proj_parm.czech; + xy_y *= this->m_proj_parm.czech; + + rho = sqrt(xy_x * xy_x + xy_y * xy_y); eps = atan2(xy_y, xy_x); - d = eps / sin(s0); - s = 2. * (atan( pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45); - u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d)); - deltav = asin(cos(s) * sin(d) / cos(u)); + d = eps / sin(S0); + 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); - lp_lon = this->m_par.lam0 - deltav / alfa; + u = asin(cos(this->m_proj_parm.ad) * sin(s) - sin(this->m_proj_parm.ad) * cos(s) * cos(d)); + deltav = asin(cos(s) * sin(d) / cos(u)); - /* ITERATION FOR lp_lat */ - fi1 = u; + lp_lon = this->m_par.lam0 - deltav / this->m_proj_parm.alpha; - ok = 0; - do - { - lp_lat = 2. * ( atan( pow( k, -1. / alfa) * - pow( tan(u / 2. + s45) , 1. / alfa) * - pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.) - ) - s45); + /* ITERATION FOR lp_lat */ + fi1 = u; - if (fabs(fi1 - lp_lat) < 0.000000000000001) ok=1; - fi1 = lp_lat; + for (i = max_iter; i ; --i) { + lp_lat = T(2) * ( atan( math::pow( this->m_proj_parm.k, T(-1) / this->m_proj_parm.alpha) * + math::pow( tan(u / T(2) + S45) , T(1) / this->m_proj_parm.alpha) * + math::pow( (T(1) + this->m_par.e * sin(fi1)) / (T(1) - this->m_par.e * sin(fi1)) , this->m_par.e / T(2)) + ) - S45); - } - while (ok==0); + if (fabs(fi1 - lp_lat) < epsilon) + break; + fi1 = lp_lat; + } + if( i == 0 ) + BOOST_THROW_EXCEPTION( projection_exception(error_non_convergent) ); lp_lon -= this->m_par.lam0; } @@ -268,32 +197,39 @@ namespace projections template <typename Parameters, typename T> inline void setup_krovak(Parameters& par, par_krovak<T>& proj_parm) { - T ts; - /* read some Parameters, - * here Latitude Truescale */ - - ts = pj_param(par.params, "rlat_ts").f; - proj_parm.C_x = ts; + T u0, n0, g; /* we want Bessel as fixed ellipsoid */ par.a = 6377397.155; par.e = sqrt(par.es = 0.006674372230614); - /* if latitude of projection center is not set, use 49d30'N */ - if (!pj_param(par.params, "tlat_0").i) - par.phi0 = 0.863937979737193; - - /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */ - /* that will correspond to using longitudes relative to greenwich */ - /* as input and output, instead of lat/long relative to Ferro */ - if (!pj_param(par.params, "tlon_0").i) - par.lam0 = 0.7417649320975901 - 0.308341501185665; - - /* if scale not set default to 0.9999 */ - if (!pj_param(par.params, "tk").i) - par.k0 = 0.9999; - - /* always the same */ + /* if latitude of projection center is not set, use 49d30'N */ + if (!pj_param_exists(par.params, "lat_0")) + par.phi0 = 0.863937979737193; + + /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */ + /* that will correspond to using longitudes relative to greenwich */ + /* as input and output, instead of lat/long relative to Ferro */ + if (!pj_param_exists(par.params, "lon_0")) + par.lam0 = 0.7417649320975901 - 0.308341501185665; + + /* if scale not set default to 0.9999 */ + if (!pj_param_exists(par.params, "k")) + par.k0 = 0.9999; + + proj_parm.czech = 1; + if( !pj_param_exists(par.params, "czech") ) + proj_parm.czech = -1; + + /* Set up shared parameters between forward and inverse */ + proj_parm.alpha = sqrt(T(1) + (par.es * math::pow(cos(par.phi0), 4)) / (T(1) - par.es)); + u0 = asin(sin(par.phi0) / proj_parm.alpha); + g = math::pow( (T(1) + par.e * sin(par.phi0)) / (T(1) - par.e * sin(par.phi0)) , proj_parm.alpha * par.e / T(2) ); + proj_parm.k = tan( u0 / 2. + S45) / math::pow(tan(par.phi0 / T(2) + S45) , proj_parm.alpha) * g; + n0 = sqrt(T(1) - par.es) / (T(1) - par.es * math::pow(sin(par.phi0), 2)); + proj_parm.n = sin(S0); + proj_parm.rho0 = par.k0 * n0 / tan(S0); + proj_parm.ad = S90 - UQ; } }} // namespace detail::krovak @@ -316,10 +252,10 @@ namespace projections \par Example \image html ex_krovak.gif */ - template <typename CalculationType, typename Parameters> - struct krovak_ellipsoid : public detail::krovak::base_krovak_ellipsoid<CalculationType, Parameters> + template <typename T, typename Parameters> + struct krovak_ellipsoid : public detail::krovak::base_krovak_ellipsoid<T, Parameters> { - inline krovak_ellipsoid(const Parameters& par) : detail::krovak::base_krovak_ellipsoid<CalculationType, Parameters>(par) + inline krovak_ellipsoid(const Parameters& par) : detail::krovak::base_krovak_ellipsoid<T, Parameters>(par) { detail::krovak::setup_krovak(this->m_par, this->m_proj_parm); } @@ -333,20 +269,20 @@ namespace projections BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::krovak, krovak_ellipsoid, krovak_ellipsoid) // Factory entry(s) - template <typename CalculationType, typename Parameters> - class krovak_entry : public detail::factory_entry<CalculationType, Parameters> + template <typename T, typename Parameters> + class krovak_entry : public detail::factory_entry<T, Parameters> { public : - virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + virtual base_v<T, Parameters>* create_new(const Parameters& par) const { - return new base_v_fi<krovak_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + return new base_v_fi<krovak_ellipsoid<T, Parameters>, T, Parameters>(par); } }; - template <typename CalculationType, typename Parameters> - inline void krovak_init(detail::base_factory<CalculationType, Parameters>& factory) + template <typename T, typename Parameters> + inline void krovak_init(detail::base_factory<T, Parameters>& factory) { - factory.add_to_factory("krovak", new krovak_entry<CalculationType, Parameters>); + factory.add_to_factory("krovak", new krovak_entry<T, Parameters>); } } // namespace detail |