diff options
Diffstat (limited to 'boost/geometry/srs/projections/proj/etmerc.hpp')
| -rw-r--r-- | boost/geometry/srs/projections/proj/etmerc.hpp | 212 |
1 files changed, 145 insertions, 67 deletions
diff --git a/boost/geometry/srs/projections/proj/etmerc.hpp b/boost/geometry/srs/projections/proj/etmerc.hpp index 479c128f11..050aafacca 100644 --- a/boost/geometry/srs/projections/proj/etmerc.hpp +++ b/boost/geometry/srs/projections/proj/etmerc.hpp @@ -1,13 +1,9 @@ -#ifndef BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP -#define BOOST_GEOMETRY_PROJECTIONS_ETMERC_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,7 +15,7 @@ // 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: @@ -43,12 +39,28 @@ // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER // DEALINGS IN THE SOFTWARE. +/* The code in this file is largly based upon procedures: + * + * Written by: Knud Poder and Karsten Engsager + * + * Based on math from: R.Koenig and K.H. Weise, "Mathematische + * Grundlagen der hoeheren Geodaesie und Kartographie, + * Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951. + * + * Modified and used here by permission of Reference Networks + * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark +*/ + +#ifndef BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP +#define BOOST_GEOMETRY_PROJECTIONS_ETMERC_HPP + #include <boost/math/special_functions/hypot.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> #include <boost/geometry/srs/projections/impl/factory_entry.hpp> +#include <boost/geometry/srs/projections/impl/function_overloads.hpp> namespace boost { namespace geometry { @@ -56,6 +68,7 @@ namespace boost { namespace geometry namespace srs { namespace par4 { struct etmerc {}; + struct utm {}; }} //namespace srs::par4 @@ -78,25 +91,8 @@ namespace projections T gtu[6]; /* Constants for geo -> transv. merc. */ }; - /* The code in this file is largly based upon procedures: - * - * Written by: Knud Poder and Karsten Engsager - * - * Based on math from: R.Koenig and K.H. Weise, "Mathematische - * Grundlagen der hoeheren Geodaesie und Kartographie, - * Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951. - * - * Modified and used here by permission of Reference Networks - * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark - */ - - - - - template <typename T> - inline T - log1py(T const& x) { /* Compute log(1+x) accurately */ + inline T log1py(T const& x) { /* Compute log(1+x) accurately */ volatile T y = 1 + x, z = y - 1; @@ -108,16 +104,14 @@ namespace projections } template <typename T> - inline T - asinhy(T const& x) { /* Compute asinh(x) accurately */ + inline T asinhy(T const& x) { /* Compute asinh(x) accurately */ T y = fabs(x); /* Enforce odd parity */ y = log1py(y * (1 + y/(boost::math::hypot(1.0, y) + 1))); return x < 0 ? -y : y; } template <typename T> - inline T - gatg(const T *p1, int len_p1, T const& B) { + inline T gatg(const T *p1, int len_p1, T const& B) { const T *p; T h = 0, h1, h2 = 0, cos_2B; @@ -127,9 +121,9 @@ namespace projections return (B + h*sin(2*B)); } + /* Complex Clenshaw summation */ template <typename T> - inline T - clenS(const T *a, int size, T const& arg_r, T const& arg_i, T *R, T *I) { + inline T clenS(const T *a, int size, T const& arg_r, T const& arg_i, T *R, T *I) { T r, i, hr, hr1, hr2, hi, hi1, hi2; T sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i; @@ -157,14 +151,15 @@ namespace projections return(*R); } + /* Real Clenshaw summation */ template <typename T> - inline T - clens(const T *a, int size, T const& arg_r) { + inline T clens(const T *a, int size, T const& arg_r) { T r, hr, hr1, hr2, cos_arg_r; const T* p = a + size; cos_arg_r = cos(arg_r); r = 2*cos_arg_r; + /* summation loop */ for (hr1 = 0, hr = *--p; a - p;) { hr2 = hr1; @@ -175,26 +170,22 @@ namespace projections } // template class, using CRTP to implement forward/inverse - template <typename CalculationType, typename Parameters> - struct base_etmerc_ellipsoid : public base_t_fi<base_etmerc_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters> + template <typename T, typename Parameters> + struct base_etmerc_ellipsoid + : public base_t_fi<base_etmerc_ellipsoid<T, Parameters>, T, Parameters> { - - typedef CalculationType geographic_type; - typedef CalculationType cartesian_type; - - par_etmerc<CalculationType> m_proj_parm; + par_etmerc<T> m_proj_parm; inline base_etmerc_ellipsoid(const Parameters& par) - : base_t_fi<base_etmerc_ellipsoid<CalculationType, Parameters>, - CalculationType, Parameters>(*this, par) {} + : base_t_fi<base_etmerc_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 { - CalculationType sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; - CalculationType Cn = lp_lat, Ce = lp_lon; + T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; + T Cn = lp_lat, Ce = lp_lon; /* ell. LAT, LNG -> Gaussian LAT, LNG */ Cn = gatg(this->m_proj_parm.cbg, PROJ_ETMERC_ORDER, Cn); @@ -206,8 +197,9 @@ namespace projections Cn = atan2(sin_Cn, cos_Ce*cos_Cn); Ce = atan2(sin_Ce*cos_Cn, boost::math::hypot(sin_Cn, cos_Cn*cos_Ce)); + /* compl. sph. N, E -> ell. norm. N, E */ - Ce = asinhy(tan(Ce)); /* Replaces: Ce = log(tan(FORTPI + Ce*0.5)); */ + Ce = asinhy(tan(Ce)); /* Replaces: Ce = log(tan(fourth_pi + Ce*0.5)); */ Cn += clenS(this->m_proj_parm.gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe); Ce += dCe; if (fabs(Ce) <= 2.623395162778) { @@ -219,19 +211,20 @@ namespace projections // 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 { - CalculationType sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; - CalculationType Cn = xy_y, Ce = xy_x; + T sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; + T Cn = xy_y, Ce = xy_x; /* normalize N, E */ Cn = (Cn - this->m_proj_parm.Zb)/this->m_proj_parm.Qn; Ce = Ce/this->m_proj_parm.Qn; + if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */ - /* norm. N, E -> compl. sph. LAT, LNG */ + /* norm. N, E -> compl. sph. LAT, LNG */ Cn += clenS(this->m_proj_parm.utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe); Ce += dCe; - Ce = atan(sinh(Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI); */ + Ce = atan(sinh(Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - fourth_pi); */ /* compl. sph. LAT -> Gaussian LAT, LNG */ sin_Cn = sin(Cn); cos_Cn = cos(Cn); @@ -254,15 +247,17 @@ namespace projections }; - // Extended Transverse Mercator template <typename Parameters, typename T> - inline void setup_etmerc(Parameters& par, par_etmerc<T>& proj_parm) + inline void setup(Parameters& par, par_etmerc<T>& proj_parm) { T f, n, np, Z; - if (par.es <= 0) - BOOST_THROW_EXCEPTION( projection_exception(-34) ); + if (par.es <= 0) { + BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) ); + } + f = par.es / (1 + sqrt(1 - par.es)); /* Replaces: f = 1 - sqrt(1-par.es); */ + /* third flattening */ np = n = f/(2 - f); @@ -270,6 +265,7 @@ namespace projections /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */ /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */ /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */ + proj_parm.cgb[0] = n*( 2 + n*(-2/3.0 + n*(-2 + n*(116/45.0 + n*(26/45.0 + n*(-2854/675.0 )))))); proj_parm.cbg[0] = n*(-2 + n*( 2/3.0 + n*( 4/3.0 + n*(-82/45.0 + n*(32/45.0 + @@ -326,13 +322,59 @@ namespace projections np *= n; proj_parm.utg[5] = np*(-20648693/638668800.0); proj_parm.gtu[5] = np*(212378941/319334400.0); + /* Gaussian latitude value of the origin latitude */ Z = gatg(proj_parm.cbg, PROJ_ETMERC_ORDER, par.phi0); + /* Origin northing minus true northing at the origin latitude */ /* i.e. true northing = N - proj_parm.Zb */ proj_parm.Zb = - proj_parm.Qn*(Z + clens(proj_parm.gtu, PROJ_ETMERC_ORDER, 2*Z)); } + // Extended Transverse Mercator + template <typename Parameters, typename T> + inline void setup_etmerc(Parameters& par, par_etmerc<T>& proj_parm) + { + setup(par, proj_parm); + } + + // Universal Transverse Mercator (UTM) + template <typename Parameters, typename T> + inline void setup_utm(Parameters& par, par_etmerc<T>& proj_parm) + { + static const T pi = detail::pi<T>(); + + int zone; + + if (par.es == 0.0) { + BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) ); + } + + par.y0 = pj_get_param_b(par.params, "south") ? 10000000. : 0.; + par.x0 = 500000.; + if (pj_param_i(par.params, "zone", zone)) /* zone input ? */ + { + if (zone > 0 && zone <= 60) + --zone; + else { + BOOST_THROW_EXCEPTION( projection_exception(error_invalid_utm_zone) ); + } + } + else /* nearest central meridian input */ + { + zone = int_floor((adjlon(par.lam0) + pi) * 30. / pi); + if (zone < 0) + zone = 0; + else if (zone >= 60) + zone = 59; + } + par.lam0 = (zone + .5) * pi / 30. - pi; + par.k0 = 0.9996; + par.phi0 = 0.; + + setup(par, proj_parm); + } + }} // namespace detail::etmerc #endif // doxygen @@ -351,37 +393,73 @@ namespace projections \par Example \image html ex_etmerc.gif */ - template <typename CalculationType, typename Parameters> - struct etmerc_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<CalculationType, Parameters> + template <typename T, typename Parameters> + struct etmerc_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters> { - inline etmerc_ellipsoid(const Parameters& par) : detail::etmerc::base_etmerc_ellipsoid<CalculationType, Parameters>(par) + inline etmerc_ellipsoid(const Parameters& par) : detail::etmerc::base_etmerc_ellipsoid<T, Parameters>(par) { detail::etmerc::setup_etmerc(this->m_par, this->m_proj_parm); } }; + /*! + \brief Universal Transverse Mercator (UTM) projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Cylindrical + - Spheroid + \par Projection parameters + - zone: UTM Zone (integer) + - south: Denotes southern hemisphere UTM zone (boolean) + \par Example + \image html ex_utm.gif + */ + template <typename T, typename Parameters> + struct utm_ellipsoid : public detail::etmerc::base_etmerc_ellipsoid<T, Parameters> + { + inline utm_ellipsoid(const Parameters& par) : detail::etmerc::base_etmerc_ellipsoid<T, Parameters>(par) + { + detail::etmerc::setup_utm(this->m_par, this->m_proj_parm); + } + }; + #ifndef DOXYGEN_NO_DETAIL namespace detail { // Static projection BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::etmerc, etmerc_ellipsoid, etmerc_ellipsoid) + BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::utm, utm_ellipsoid, utm_ellipsoid) // Factory entry(s) - template <typename CalculationType, typename Parameters> - class etmerc_entry : public detail::factory_entry<CalculationType, Parameters> + template <typename T, typename Parameters> + class etmerc_entry : public detail::factory_entry<T, Parameters> + { + public : + virtual base_v<T, Parameters>* create_new(const Parameters& par) const + { + return new base_v_fi<etmerc_ellipsoid<T, Parameters>, T, Parameters>(par); + } + }; + + template <typename T, typename Parameters> + class utm_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<etmerc_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + return new base_v_fi<utm_ellipsoid<T, Parameters>, T, Parameters>(par); } }; - template <typename CalculationType, typename Parameters> - inline void etmerc_init(detail::base_factory<CalculationType, Parameters>& factory) + template <typename T, typename Parameters> + inline void etmerc_init(detail::base_factory<T, Parameters>& factory) { - factory.add_to_factory("etmerc", new etmerc_entry<CalculationType, Parameters>); + factory.add_to_factory("etmerc", new etmerc_entry<T, Parameters>); + factory.add_to_factory("utm", new utm_entry<T, Parameters>); } } // namespace detail |