projections/proj/cass.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: 9.0.0
0020
0021 // Original copyright notice:
0022
0023 // Permission is hereby granted, free of charge, to any person obtaining a
0024 // copy of this software and associated documentation files (the "Software"),
0025 // to deal in the Software without restriction, including without limitation
0026 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0027 // and/or sell copies of the Software, and to permit persons to whom the
0028 // Software is furnished to do so, subject to the following conditions:
0029
0030 // The above copyright notice and this permission notice shall be included
0031 // in all copies or substantial portions of the Software.
0032
0033 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0034 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0035 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0036 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0037 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0038 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0039 // DEALINGS IN THE SOFTWARE.
0040
0041 #ifndef BOOST_GEOMETRY_PROJECTIONS_CASS_HPP
0042 #define BOOST_GEOMETRY_PROJECTIONS_CASS_HPP
0043
0044 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0045 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0046 #include <boost/geometry/srs/projections/impl/projects.hpp>
0047 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0048 #include <boost/geometry/srs/projections/impl/pj_generic_inverse.hpp>
0049 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
0050 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0051
0052 namespace boost { namespace geometry
0053 {
0054
0055 namespace projections
0056 {
0057     #ifndef DOXYGEN_NO_DETAIL
0058     namespace detail { namespace cass
0059     {
0060
0061             //static const double epsilon10 = 1e-10;
0062             //static const double C1 = .16666666666666666666;
0063             //static const double C2 = .00833333333333333333;
0064             //static const double C3 = .04166666666666666666;
0065             //static const double C4 = .33333333333333333333;
0066             //static const double C5 = .06666666666666666666;
0067
0068             template <typename T>
0069             inline T C1() { return .16666666666666666666666666666666666666; }
0070             template <typename T>
0071             inline T C2() { return .00833333333333333333333333333333333333; }
0072             template <typename T>
0073             inline T C3() { return .04166666666666666666666666666666666666; }
0074             template <typename T>
0075             inline T C4() { return .33333333333333333333333333333333333333; }
0076             template <typename T>
0077             inline T C5() { return .06666666666666666666666666666666666666; }
0078
0079             template <typename T>
0080             struct par_cass
0081             {
0082                 T m0;
0083                 detail::en<T> en;
0084                 bool hyperbolic;
0085             };
0086
0087             template <typename T, typename Parameters>
0088             struct base_cass_ellipsoid
0089             {
0090                 par_cass<T> m_proj_parm;
0091
0092                 // FORWARD(e_forward)  ellipsoid
0093                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0094                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat,
0095                                 T& xy_x, T& xy_y) const
0096                 {
0097                     static T const C1 = cass::C1<T>();
0098                     static T const C2 = cass::C2<T>();
0099                     static T const C3 = cass::C3<T>();
0100
0101                     T sinlplat = sin(lp_lat);
0102                     T c = cos(lp_lat);
0103                     xy_y = pj_mlfn(lp_lat, sinlplat, c, this->m_proj_parm.en);
0104                     T n_square = 1./(1. - par.es * sinlplat * sinlplat);
0105                     T n = sqrt(n_square);
0106                     T tn = tan(lp_lat);
0107                     T t = tn * tn;
0108                     T a1 = lp_lon * c;
0109                     c *= par.es * c / (1 - par.es);
0110                     T a2 = a1 * a1;
0111                     xy_x = n * a1 * (1. - a2 * t *
0112                         (C1 - (8. - t + 8. * c) * a2 * C2));
0113                     xy_y -= this->m_proj_parm.m0 - n * tn * a2 *
0114                         (.5 + (5. - t + 6. * c) * a2 * C3);
0115                     if ( this->m_proj_parm.hyperbolic )
0116                     {
0117                         T const rho = n_square * (1. - par.es) * n;
0118                         xy_y -= xy_y * xy_y * xy_y / (6 * rho * n);
0119                     }
0120                 }
0121
0122                 // INVERSE(e_inverse)  ellipsoid
0123                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0124                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y,
0125                                 T& lp_lon, T& lp_lat) const
0126                 {
0127                     static const T C3 = cass::C3<T>();
0128                     static const T C4 = cass::C4<T>();
0129                     static const T C5 = cass::C5<T>();
0130
0131                     T ph1;
0132
0133                     ph1 = pj_inv_mlfn(this->m_proj_parm.m0 + xy_y, par.es, this->m_proj_parm.en);
0134                     T tn = tan(ph1); T t = tn * tn;
0135                     T n = sin(ph1);
0136                     T r = 1. / (1. - par.es * n * n);
0137                     n = sqrt(r);
0138                     r *= (1. - par.es) * n;
0139                     T dd = xy_x / n;
0140                     T d2 = dd * dd;
0141                     lp_lat = ph1 - (n * tn / r) * d2 *
0142                         (.5 - (1. + 3. * t) * d2 * C3);
0143                     lp_lon = dd * (1. + t * d2 *
0144                         (-C4 + (1. + 3. * t) * d2 * C5)) / cos(ph1);
0145                     if ( this->m_proj_parm.hyperbolic )
0146                     {
0147                         // EPSG guidance note 7-2 suggests a custom approximation for the
0148                         // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the
0149                         // generic inversion method
0150                         pj_generic_inverse_2d(xy_x, xy_y, par, this, lp_lat, lp_lon);
0151                     }
0152                 }
0153
0154                 static inline std::string get_name()
0155                 {
0156                     return "cass_ellipsoid";
0157                 }
0158
0159             };
0160
0161             template <typename T, typename Parameters>
0162             struct base_cass_spheroid
0163             {
0164                 par_cass<T> m_proj_parm;
0165
0166                 // FORWARD(s_forward)  spheroid
0167                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0168                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat,
0169                                 T& xy_x, T& xy_y) const
0170                 {
0171                     xy_x = asin(cos(lp_lat) * sin(lp_lon));
0172                     xy_y = atan2(tan(lp_lat) , cos(lp_lon)) - par.phi0;
0173                 }
0174
0175                 // INVERSE(s_inverse)  spheroid
0176                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0177                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y,
0178                                 T& lp_lon, T& lp_lat) const
0179                 {
0180                     T dd = xy_y + par.phi0;
0181                     lp_lat = asin(sin(dd) * cos(xy_x));
0182                     lp_lon = atan2(tan(xy_x), cos(dd));
0183                 }
0184
0185                 static inline std::string get_name()
0186                 {
0187                     return "cass_spheroid";
0188                 }
0189
0190             };
0191
0192             // Cassini
0193             template <typename Params, typename Parameters, typename T>
0194             inline void setup_cass(Params const& params, Parameters& par, par_cass<T>& proj_parm)
0195             {
0196                 if (par.es) {
0197                     proj_parm.en = pj_enfn<T>(par.es);
0198                     proj_parm.m0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en);
0199                 } else {
0200                 }
0201                 proj_parm.hyperbolic = false;
0202                 if (pj_param_exists<srs::spar::hyperbolic>(params, "hyperbolic", srs::dpar::hyperbolic))
0203                     proj_parm.hyperbolic = true;
0204
0205             }
0206
0207     }} // namespace detail::cass
0208     #endif // doxygen
0209
0210     /*!
0211         \brief Cassini projection
0212         \ingroup projections
0213         \tparam Geographic latlong point type
0214         \tparam Cartesian xy point type
0215         \tparam Parameters parameter type
0216         \par Projection characteristics
0217          - Cylindrical
0218          - Spheroid
0219          - Ellipsoid
0220         \par Example
0221         \image html ex_cass.gif
0222     */
0223     template <typename T, typename Parameters>
0224     struct cass_ellipsoid : public detail::cass::base_cass_ellipsoid<T, Parameters>
0225     {
0226         template <typename Params>
0227         inline cass_ellipsoid(Params const& params, Parameters const& par)
0228         {
0229             detail::cass::setup_cass(params, par, this->m_proj_parm);
0230         }
0231     };
0232
0233     /*!
0234         \brief Cassini projection
0235         \ingroup projections
0236         \tparam Geographic latlong point type
0237         \tparam Cartesian xy point type
0238         \tparam Parameters parameter type
0239         \par Projection characteristics
0240          - Cylindrical
0241          - Spheroid
0242          - Ellipsoid
0243         \par Example
0244         \image html ex_cass.gif
0245     */
0246     template <typename T, typename Parameters>
0247     struct cass_spheroid : public detail::cass::base_cass_spheroid<T, Parameters>
0248     {
0249         template <typename Params>
0250         inline cass_spheroid(Params const& params, Parameters const& par)
0251         {
0252             detail::cass::setup_cass(params, par, this->m_proj_parm);
0253         }
0254     };
0255
0256     #ifndef DOXYGEN_NO_DETAIL
0257     namespace detail
0258     {
0259
0260         // Static projection
0261         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_cass, cass_spheroid,
0262                                                                 cass_ellipsoid)
0263
0264         // Factory entry(s)
0265         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(cass_entry, cass_spheroid, cass_ellipsoid)
0266
0267         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(cass_init)
0268         {
0269             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(cass, cass_entry);
0270         }
0271
0272     } // namespace detail
0273     #endif // doxygen
0274
0275 } // namespace projections
0276
0277 }} // namespace boost::geometry
0278
0279 #endif // BOOST_GEOMETRY_PROJECTIONS_CASS_HPP
0280