boost/geometry/srs/projections/proj/tmerc.hpp
#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