diff options
Diffstat (limited to 'boost/geometry/srs/projections/proj/tmerc.hpp')
| -rw-r--r-- | boost/geometry/srs/projections/proj/tmerc.hpp | 495 |
1 files changed, 495 insertions, 0 deletions
diff --git a/boost/geometry/srs/projections/proj/tmerc.hpp b/boost/geometry/srs/projections/proj/tmerc.hpp new file mode 100644 index 0000000000..6ed2902142 --- /dev/null +++ b/boost/geometry/srs/projections/proj/tmerc.hpp @@ -0,0 +1,495 @@ +#ifndef BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP +#define BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP + +// Boost.Geometry - extensions-gis-projections (based on PROJ4) +// This file is automatically generated. DO NOT EDIT. + +// 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. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +// This file is converted from PROJ4, http://trac.osgeo.org/proj +// PROJ4 is originally written by Gerald Evenden (then of the USGS) +// PROJ4 is maintained by Frank Warmerdam +// PROJ4 is converted to Boost.Geometry by Barend Gehrels + +// Last updated version of proj: 4.9.1 + +// Original copyright notice: + +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: + +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. + +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. + +#include <boost/geometry/util/math.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> +#include <boost/geometry/srs/projections/impl/pj_mlfn.hpp> + + +namespace boost { namespace geometry +{ + +namespace srs { namespace par4 +{ + struct tmerc {}; + struct utm {}; + +}} //namespace srs::par4 + +namespace projections +{ + #ifndef DOXYGEN_NO_DETAIL + namespace detail { namespace tmerc + { + + static const double EPS10 = 1.e-10; + //static const double FC1 = 1.; + //static const double FC2 = .5; + //static const double FC3 = .16666666666666666666; + //static const double FC4 = .08333333333333333333; + //static const double FC5 = .05; + //static const double FC6 = .03333333333333333333; + //static const double FC7 = .02380952380952380952; + //static const double FC8 = .01785714285714285714; + + template <typename T> + inline T FC1() { return 1.; } + template <typename T> + inline T FC2() { return .5; } + template <typename T> + inline T FC3() { return .16666666666666666666666666666666666666; } + template <typename T> + inline T FC4() { return .08333333333333333333333333333333333333; } + template <typename T> + inline T FC5() { return .05; } + template <typename T> + inline T FC6() { return .03333333333333333333333333333333333333; } + template <typename T> + inline T FC7() { return .02380952380952380952380952380952380952; } + template <typename T> + inline T FC8() { return .01785714285714285714285714285714285714; } + + template <typename T> + struct par_tmerc + { + T esp; + T ml0; + T en[EN_SIZE]; + }; + + // template class, using CRTP to implement forward/inverse + template <typename CalculationType, typename Parameters> + struct base_tmerc_ellipsoid : public base_t_fi<base_tmerc_ellipsoid<CalculationType, Parameters>, + CalculationType, Parameters> + { + + typedef CalculationType geographic_type; + typedef CalculationType cartesian_type; + + par_tmerc<CalculationType> m_proj_parm; + + inline base_tmerc_ellipsoid(const Parameters& par) + : base_t_fi<base_tmerc_ellipsoid<CalculationType, Parameters>, + CalculationType, Parameters>(*this, par) {} + + // FORWARD(e_forward) ellipse + // 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 + { + static const CalculationType HALFPI = detail::HALFPI<CalculationType>(); + static const CalculationType FC1 = tmerc::FC1<CalculationType>(); + static const CalculationType FC2 = tmerc::FC2<CalculationType>(); + static const CalculationType FC3 = tmerc::FC3<CalculationType>(); + static const CalculationType FC4 = tmerc::FC4<CalculationType>(); + static const CalculationType FC5 = tmerc::FC5<CalculationType>(); + static const CalculationType FC6 = tmerc::FC6<CalculationType>(); + static const CalculationType FC7 = tmerc::FC7<CalculationType>(); + static const CalculationType FC8 = tmerc::FC8<CalculationType>(); + + CalculationType al, als, n, cosphi, sinphi, t; + + /* + * Fail if our longitude is more than 90 degrees from the + * central meridian since the results are essentially garbage. + * Is error -20 really an appropriate return value? + * + * http://trac.osgeo.org/proj/ticket/5 + */ + if( lp_lon < -HALFPI || lp_lon > HALFPI ) + { + xy_x = HUGE_VAL; + xy_y = HUGE_VAL; + BOOST_THROW_EXCEPTION( projection_exception(-14) ); + return; + } + + sinphi = sin(lp_lat); + cosphi = cos(lp_lat); + t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; + t *= t; + al = cosphi * lp_lon; + als = al * al; + al /= sqrt(1. - this->m_par.es * sinphi * sinphi); + n = this->m_proj_parm.esp * cosphi * cosphi; + xy_x = this->m_par.k0 * al * (FC1 + + FC3 * als * (1. - t + n + + FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t) + + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) ) + ))); + xy_y = this->m_par.k0 * (pj_mlfn(lp_lat, sinphi, cosphi, this->m_proj_parm.en) - this->m_proj_parm.ml0 + + sinphi * al * lp_lon * FC2 * ( 1. + + FC4 * als * (5. - t + n * (9. + 4. * n) + + FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t) + + FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) ) + )))); + } + + // 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 + { + static const CalculationType HALFPI = detail::HALFPI<CalculationType>(); + static const CalculationType FC1 = tmerc::FC1<CalculationType>(); + static const CalculationType FC2 = tmerc::FC2<CalculationType>(); + static const CalculationType FC3 = tmerc::FC3<CalculationType>(); + static const CalculationType FC4 = tmerc::FC4<CalculationType>(); + static const CalculationType FC5 = tmerc::FC5<CalculationType>(); + static const CalculationType FC6 = tmerc::FC6<CalculationType>(); + static const CalculationType FC7 = tmerc::FC7<CalculationType>(); + static const CalculationType FC8 = tmerc::FC8<CalculationType>(); + + CalculationType n, con, cosphi, d, ds, sinphi, t; + + lp_lat = pj_inv_mlfn(this->m_proj_parm.ml0 + xy_y / this->m_par.k0, this->m_par.es, this->m_proj_parm.en); + if (fabs(lp_lat) >= HALFPI) { + lp_lat = xy_y < 0. ? -HALFPI : HALFPI; + lp_lon = 0.; + } else { + sinphi = sin(lp_lat); + cosphi = cos(lp_lat); + t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.; + n = this->m_proj_parm.esp * cosphi * cosphi; + d = xy_x * sqrt(con = 1. - this->m_par.es * sinphi * sinphi) / this->m_par.k0; + con *= t; + t *= t; + ds = d * d; + lp_lat -= (con * ds / (1.-this->m_par.es)) * FC2 * (1. - + ds * FC4 * (5. + t * (3. - 9. * n) + n * (1. - 4 * n) - + ds * FC6 * (61. + t * (90. - 252. * n + + 45. * t) + 46. * n + - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) ) + ))); + lp_lon = d*(FC1 - + ds*FC3*( 1. + 2.*t + n - + ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n + - ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) ) + ))) / cosphi; + } + } + + static inline std::string get_name() + { + return "tmerc_ellipsoid"; + } + + }; + + // template class, using CRTP to implement forward/inverse + template <typename CalculationType, typename Parameters> + struct base_tmerc_spheroid : public base_t_fi<base_tmerc_spheroid<CalculationType, Parameters>, + CalculationType, Parameters> + { + + typedef CalculationType geographic_type; + typedef CalculationType cartesian_type; + + par_tmerc<CalculationType> m_proj_parm; + + inline base_tmerc_spheroid(const Parameters& par) + : base_t_fi<base_tmerc_spheroid<CalculationType, Parameters>, + CalculationType, Parameters>(*this, par) {} + + // FORWARD(s_forward) sphere + // 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 + { + static const CalculationType HALFPI = detail::HALFPI<CalculationType>(); + + CalculationType b, cosphi; + + /* + * Fail if our longitude is more than 90 degrees from the + * central meridian since the results are essentially garbage. + * Is error -20 really an appropriate return value? + * + * http://trac.osgeo.org/proj/ticket/5 + */ + if( lp_lon < -HALFPI || lp_lon > HALFPI ) + { + xy_x = HUGE_VAL; + xy_y = HUGE_VAL; + BOOST_THROW_EXCEPTION( projection_exception(-14) ); + return; + } + + cosphi = cos(lp_lat); + b = cosphi * sin(lp_lon); + if (fabs(fabs(b) - 1.) <= EPS10) + BOOST_THROW_EXCEPTION( projection_exception(-20) ); + + xy_x = this->m_proj_parm.ml0 * log((1. + b) / (1. - b)); + xy_y = cosphi * cos(lp_lon) / sqrt(1. - b * b); + + b = fabs( xy_y ); + if (b >= 1.) { + if ((b - 1.) > EPS10) + BOOST_THROW_EXCEPTION( projection_exception(-20) ); + else xy_y = 0.; + } else + xy_y = acos(xy_y); + + if (lp_lat < 0.) + xy_y = -xy_y; + xy_y = this->m_proj_parm.esp * (xy_y - this->m_par.phi0); + } + + // INVERSE(s_inverse) sphere + // 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 + { + CalculationType h, g; + + h = exp(xy_x / this->m_proj_parm.esp); + g = .5 * (h - 1. / h); + h = cos(this->m_par.phi0 + xy_y / this->m_proj_parm.esp); + lp_lat = asin(sqrt((1. - h * h) / (1. + g * g))); + if (xy_y < 0.) lp_lat = -lp_lat; + lp_lon = (g || h) ? atan2(g, h) : 0.; + } + + static inline std::string get_name() + { + return "tmerc_spheroid"; + } + + }; + + template <typename Parameters, typename T> + inline void setup(Parameters& par, par_tmerc<T>& proj_parm) /* general initialization */ + { + if (par.es) { + if (!pj_enfn(par.es, proj_parm.en)) + BOOST_THROW_EXCEPTION( projection_exception(0) ); + proj_parm.ml0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en); + proj_parm.esp = par.es / (1. - par.es); + } else { + proj_parm.esp = par.k0; + proj_parm.ml0 = .5 * proj_parm.esp; + } + } + + + // Transverse Mercator + template <typename Parameters, typename T> + inline void setup_tmerc(Parameters& par, par_tmerc<T>& proj_parm) + { + setup(par, proj_parm); + } + + // Universal Transverse Mercator (UTM) + template <typename Parameters, typename T> + inline void setup_utm(Parameters& par, par_tmerc<T>& proj_parm) + { + static const T ONEPI = detail::ONEPI<T>(); + + int zone; + + par.y0 = pj_param(par.params, "bsouth").i ? 10000000. : 0.; + par.x0 = 500000.; + if (pj_param(par.params, "tzone").i) /* zone input ? */ + if ((zone = pj_param(par.params, "izone").i) > 0 && zone <= 60) + --zone; + else + BOOST_THROW_EXCEPTION( projection_exception(-35) ); + else /* nearest central meridian input */ + if ((zone = int_floor((adjlon(par.lam0) + ONEPI) * 30. / ONEPI)) < 0) + zone = 0; + else if (zone >= 60) + zone = 59; + par.lam0 = (zone + .5) * ONEPI / 30. - ONEPI; + par.k0 = 0.9996; + par.phi0 = 0.; + setup(par, proj_parm); + } + + }} // namespace detail::tmerc + #endif // doxygen + + /*! + \brief Transverse Mercator projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Cylindrical + - Spheroid + - Ellipsoid + \par Example + \image html ex_tmerc.gif + */ + template <typename CalculationType, typename Parameters> + struct tmerc_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid<CalculationType, Parameters> + { + inline tmerc_ellipsoid(const Parameters& par) : detail::tmerc::base_tmerc_ellipsoid<CalculationType, Parameters>(par) + { + detail::tmerc::setup_tmerc(this->m_par, this->m_proj_parm); + } + }; + + /*! + \brief Transverse Mercator projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Cylindrical + - Spheroid + - Ellipsoid + \par Example + \image html ex_tmerc.gif + */ + template <typename CalculationType, typename Parameters> + struct tmerc_spheroid : public detail::tmerc::base_tmerc_spheroid<CalculationType, Parameters> + { + inline tmerc_spheroid(const Parameters& par) : detail::tmerc::base_tmerc_spheroid<CalculationType, Parameters>(par) + { + detail::tmerc::setup_tmerc(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 CalculationType, typename Parameters> + struct utm_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid<CalculationType, Parameters> + { + inline utm_ellipsoid(const Parameters& par) : detail::tmerc::base_tmerc_ellipsoid<CalculationType, Parameters>(par) + { + detail::tmerc::setup_utm(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 CalculationType, typename Parameters> + struct utm_spheroid : public detail::tmerc::base_tmerc_spheroid<CalculationType, Parameters> + { + inline utm_spheroid(const Parameters& par) : detail::tmerc::base_tmerc_spheroid<CalculationType, Parameters>(par) + { + detail::tmerc::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::tmerc, tmerc_spheroid, tmerc_ellipsoid) + BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::utm, utm_spheroid, utm_ellipsoid) + + // Factory entry(s) - dynamic projection + template <typename CalculationType, typename Parameters> + class tmerc_entry : public detail::factory_entry<CalculationType, Parameters> + { + public : + virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + { + if (par.es) + return new base_v_fi<tmerc_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + else + return new base_v_fi<tmerc_spheroid<CalculationType, Parameters>, CalculationType, Parameters>(par); + } + }; + + template <typename CalculationType, typename Parameters> + class utm_entry : public detail::factory_entry<CalculationType, Parameters> + { + public : + virtual base_v<CalculationType, Parameters>* create_new(const Parameters& par) const + { + if (par.es) + return new base_v_fi<utm_ellipsoid<CalculationType, Parameters>, CalculationType, Parameters>(par); + else + return new base_v_fi<utm_spheroid<CalculationType, Parameters>, CalculationType, Parameters>(par); + } + }; + + template <typename CalculationType, typename Parameters> + inline void tmerc_init(detail::base_factory<CalculationType, Parameters>& factory) + { + factory.add_to_factory("tmerc", new tmerc_entry<CalculationType, Parameters>); + factory.add_to_factory("utm", new utm_entry<CalculationType, Parameters>); + } + + } // namespace detail + #endif // doxygen + +} // namespace projections + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_PROJECTIONS_TMERC_HPP + |