projections/proj/imw_p.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.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Permission is hereby granted, free of charge, to any person obtaining a
0023 // copy of this software and associated documentation files (the "Software"),
0024 // to deal in the Software without restriction, including without limitation
0025 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0026 // and/or sell copies of the Software, and to permit persons to whom the
0027 // Software is furnished to do so, subject to the following conditions:
0028
0029 // The above copyright notice and this permission notice shall be included
0030 // in all copies or substantial portions of the Software.
0031
0032 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0033 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0034 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0035 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0036 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0037 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0038 // DEALINGS IN THE SOFTWARE.
0039
0040 #ifndef BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
0042
0043 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0044 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0045 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0046 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
0047 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0048 #include <boost/geometry/srs/projections/impl/projects.hpp>
0049
0050 #include <boost/geometry/util/math.hpp>
0051
0052 namespace boost { namespace geometry
0053 {
0054
0055 namespace projections
0056 {
0057     #ifndef DOXYGEN_NO_DETAIL
0058     namespace detail { namespace imw_p
0059     {
0060
0061             static const double tolerance = 1e-10;
0062             static const double epsilon = 1e-10;
0063
0064             template <typename T>
0065             struct point_xy { T x, y; }; // specific for IMW_P
0066
0067             enum mode_type {
0068                 none_is_zero  =  0, /* phi_1 and phi_2 != 0 */
0069                 phi_1_is_zero =  1, /* phi_1 = 0 */
0070                 phi_2_is_zero = -1  /* phi_2 = 0 */
0071             };
0072
0073             template <typename T>
0074             struct par_imw_p
0075             {
0076                 T    P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
0077                 T    phi_1, phi_2, lam_1;
0078                 detail::en<T> en;
0079                 mode_type mode;
0080             };
0081
0082             template <typename Params, typename T>
0083             inline int phi12(Params const& params,
0084                              par_imw_p<T> & proj_parm, T *del, T *sig)
0085             {
0086                 int err = 0;
0087
0088                 if (!pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi_1) ||
0089                     !pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi_2)) {
0090                     err = -41;
0091                 } else {
0092                     //proj_parm.phi_1 = pj_get_param_r(par.params, "lat_1"); // set above
0093                     //proj_parm.phi_2 = pj_get_param_r(par.params, "lat_2"); // set above
0094                     *del = 0.5 * (proj_parm.phi_2 - proj_parm.phi_1);
0095                     *sig = 0.5 * (proj_parm.phi_2 + proj_parm.phi_1);
0096                     err = (fabs(*del) < epsilon || fabs(*sig) < epsilon) ? -42 : 0;
0097                 }
0098                 return err;
0099             }
0100             template <typename Parameters, typename T>
0101             inline point_xy<T> loc_for(T const& lp_lam, T const& lp_phi,
0102                                        Parameters const& par,
0103                                        par_imw_p<T> const& proj_parm,
0104                                        T *yc)
0105             {
0106                 point_xy<T> xy;
0107
0108                 if (lp_phi == 0.0) {
0109                     xy.x = lp_lam;
0110                     xy.y = 0.;
0111                 } else {
0112                     T xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
0113
0114                     sp = sin(lp_phi);
0115                     m = pj_mlfn(lp_phi, sp, cos(lp_phi), proj_parm.en);
0116                     xa = proj_parm.Pp + proj_parm.Qp * m;
0117                     ya = proj_parm.P + proj_parm.Q * m;
0118                     R = 1. / (tan(lp_phi) * sqrt(1. - par.es * sp * sp));
0119                     C = sqrt(R * R - xa * xa);
0120                     if (lp_phi < 0.) C = - C;
0121                     C += ya - R;
0122                     if (proj_parm.mode == phi_2_is_zero) {
0123                         xb = lp_lam;
0124                         yb = proj_parm.C2;
0125                     } else {
0126                         t = lp_lam * proj_parm.sphi_2;
0127                         xb = proj_parm.R_2 * sin(t);
0128                         yb = proj_parm.C2 + proj_parm.R_2 * (1. - cos(t));
0129                     }
0130                     if (proj_parm.mode == phi_1_is_zero) {
0131                         xc = lp_lam;
0132                         *yc = 0.;
0133                     } else {
0134                         t = lp_lam * proj_parm.sphi_1;
0135                         xc = proj_parm.R_1 * sin(t);
0136                         *yc = proj_parm.R_1 * (1. - cos(t));
0137                     }
0138                     D = (xb - xc)/(yb - *yc);
0139                     B = xc + D * (C + R - *yc);
0140                     xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
0141                     if (lp_phi > 0)
0142                         xy.x = - xy.x;
0143                     xy.x = (B + xy.x) / (1. + D * D);
0144                     xy.y = sqrt(R * R - xy.x * xy.x);
0145                     if (lp_phi > 0)
0146                         xy.y = - xy.y;
0147                     xy.y += C + R;
0148                 }
0149                 return (xy);
0150             }
0151             template <typename Parameters, typename T>
0152             inline void xy(Parameters const& par, par_imw_p<T> const& proj_parm,
0153                            T const& phi,
0154                            T *x, T *y, T *sp, T *R)
0155             {
0156                 T F;
0157
0158                 *sp = sin(phi);
0159                 *R = 1./(tan(phi) * sqrt(1. - par.es * *sp * *sp ));
0160                 F = proj_parm.lam_1 * *sp;
0161                 *y = *R * (1 - cos(F));
0162                 *x = *R * sin(F);
0163             }
0164
0165             template <typename T, typename Parameters>
0166             struct base_imw_p_ellipsoid
0167             {
0168                 par_imw_p<T> m_proj_parm;
0169
0170                 // FORWARD(e_forward)  ellipsoid
0171                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0172                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0173                 {
0174                     T yc = 0;
0175                     point_xy<T> xy = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
0176                     xy_x = xy.x; xy_y = xy.y;
0177                 }
0178
0179                 // INVERSE(e_inverse)  ellipsoid
0180                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0181                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0182                 {
0183                     point_xy<T> t;
0184                     T yc = 0.0;
0185                     int i = 0;
0186                     const int n_max_iter = 1000; /* Arbitrarily choosen number... */
0187
0188                     lp_lat = this->m_proj_parm.phi_2;
0189                     lp_lon = xy_x / cos(lp_lat);
0190                     do {
0191                         t = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
0192                         lp_lat = ((lp_lat - this->m_proj_parm.phi_1) * (xy_y - yc) / (t.y - yc)) + this->m_proj_parm.phi_1;
0193                         lp_lon = lp_lon * xy_x / t.x;
0194                         i++;
0195                     } while (i < n_max_iter &&
0196                              (fabs(t.x - xy_x) > tolerance || fabs(t.y - xy_y) > tolerance));
0197
0198                     if( i == n_max_iter )
0199                     {
0200                         lp_lon = lp_lat = HUGE_VAL;
0201                     }
0202                 }
0203
0204                 static inline std::string get_name()
0205                 {
0206                     return "imw_p_ellipsoid";
0207                 }
0208
0209             };
0210
0211             // International Map of the World Polyconic
0212             template <typename Params, typename Parameters, typename T>
0213             inline void setup_imw_p(Params const& params, Parameters const& par, par_imw_p<T>& proj_parm)
0214             {
0215                 T del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
0216                 int err;
0217
0218                 proj_parm.en = pj_enfn<T>(par.es);
0219                 if( (err = phi12(params, proj_parm, &del, &sig)) != 0)
0220                     BOOST_THROW_EXCEPTION( projection_exception(err) );
0221                 if (proj_parm.phi_2 < proj_parm.phi_1) { /* make sure proj_parm.phi_1 most southerly */
0222                     del = proj_parm.phi_1;
0223                     proj_parm.phi_1 = proj_parm.phi_2;
0224                     proj_parm.phi_2 = del;
0225                 }
0226                 if (pj_param_r<srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1, proj_parm.lam_1)) {
0227                     /* empty */
0228                 } else { /* use predefined based upon latitude */
0229                     sig = fabs(sig * geometry::math::r2d<T>());
0230                     if (sig <= 60)      sig = 2.;
0231                     else if (sig <= 76) sig = 4.;
0232                     else                sig = 8.;
0233                     proj_parm.lam_1 = sig * geometry::math::d2r<T>();
0234                 }
0235                 proj_parm.mode = none_is_zero;
0236                 if (proj_parm.phi_1 != 0.0)
0237                     xy(par, proj_parm, proj_parm.phi_1, &x1, &y1, &proj_parm.sphi_1, &proj_parm.R_1);
0238                 else {
0239                     proj_parm.mode = phi_1_is_zero;
0240                     y1 = 0.;
0241                     x1 = proj_parm.lam_1;
0242                 }
0243                 if (proj_parm.phi_2 != 0.0)
0244                     xy(par, proj_parm, proj_parm.phi_2, &x2, &T2, &proj_parm.sphi_2, &proj_parm.R_2);
0245                 else {
0246                     proj_parm.mode = phi_2_is_zero;
0247                     T2 = 0.;
0248                     x2 = proj_parm.lam_1;
0249                 }
0250                 m1 = pj_mlfn(proj_parm.phi_1, proj_parm.sphi_1, cos(proj_parm.phi_1), proj_parm.en);
0251                 m2 = pj_mlfn(proj_parm.phi_2, proj_parm.sphi_2, cos(proj_parm.phi_2), proj_parm.en);
0252                 t = m2 - m1;
0253                 s = x2 - x1;
0254                 y2 = sqrt(t * t - s * s) + y1;
0255                 proj_parm.C2 = y2 - T2;
0256                 t = 1. / t;
0257                 proj_parm.P = (m2 * y1 - m1 * y2) * t;
0258                 proj_parm.Q = (y2 - y1) * t;
0259                 proj_parm.Pp = (m2 * x1 - m1 * x2) * t;
0260                 proj_parm.Qp = (x2 - x1) * t;
0261             }
0262
0263     }} // namespace detail::imw_p
0264     #endif // doxygen
0265
0266     /*!
0267         \brief International Map of the World Polyconic projection
0268         \ingroup projections
0269         \tparam Geographic latlong point type
0270         \tparam Cartesian xy point type
0271         \tparam Parameters parameter type
0272         \par Projection characteristics
0273          - Mod. Polyconic
0274          - Ellipsoid
0275         \par Projection parameters
0276          - lat_1: Latitude of first standard parallel
0277          - lat_2: Latitude of second standard parallel
0278          - lon_1 (degrees)
0279         \par Example
0280         \image html ex_imw_p.gif
0281     */
0282     template <typename T, typename Parameters>
0283     struct imw_p_ellipsoid : public detail::imw_p::base_imw_p_ellipsoid<T, Parameters>
0284     {
0285         template <typename Params>
0286         inline imw_p_ellipsoid(Params const& params, Parameters const& par)
0287         {
0288             detail::imw_p::setup_imw_p(params, par, this->m_proj_parm);
0289         }
0290     };
0291
0292     #ifndef DOXYGEN_NO_DETAIL
0293     namespace detail
0294     {
0295
0296         // Static projection
0297         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_imw_p, imw_p_ellipsoid)
0298
0299         // Factory entry(s)
0300         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(imw_p_entry, imw_p_ellipsoid)
0301
0302         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(imw_p_init)
0303         {
0304             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(imw_p, imw_p_entry)
0305         }
0306
0307     } // namespace detail
0308     #endif // doxygen
0309
0310 } // namespace projections
0311
0312 }} // namespace boost::geometry
0313
0314 #endif // BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
0315