boost/geometry/srs/projections/proj/imw_p.hpp

   1 // Boost.Geometry - gis-projections (based on PROJ4)
   2
   3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4
   5 // This file was modified by Oracle on 2017, 2018, 2019.
   6 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
   7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
   8
   9 // Use, modification and distribution is subject to the Boost Software License,
  10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  11 // http://www.boost.org/LICENSE_1_0.txt)
  12
  13 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  14 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  15 // PROJ4 is maintained by Frank Warmerdam
  16 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
  17
  18 // Last updated version of proj: 5.0.0
  19
  20 // Original copyright notice:
  21
  22 // Permission is hereby granted, free of charge, to any person obtaining a
  23 // copy of this software and associated documentation files (the "Software"),
  24 // to deal in the Software without restriction, including without limitation
  25 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  26 // and/or sell copies of the Software, and to permit persons to whom the
  27 // Software is furnished to do so, subject to the following conditions:
  28
  29 // The above copyright notice and this permission notice shall be included
  30 // in all copies or substantial portions of the Software.
  31
  32 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  33 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  34 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  35 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  36 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  37 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  38 // DEALINGS IN THE SOFTWARE.
  39
  40 #ifndef BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
  41 #define BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
  42
  43 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  44 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  45 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  46 #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
  47 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  48 #include <boost/geometry/srs/projections/impl/projects.hpp>
  49
  50 #include <boost/geometry/util/math.hpp>
  51
  52 namespace boost { namespace geometry
  53 {
  54
  55 namespace projections
  56 {
  57     #ifndef DOXYGEN_NO_DETAIL
  58     namespace detail { namespace imw_p
  59     {
  60
  61             static const double tolerance = 1e-10;
  62             static const double epsilon = 1e-10;
  63
  64             template <typename T>
  65             struct point_xy { T x, y; }; // specific for IMW_P
  66
  67             enum mode_type {
  68                 none_is_zero  =  0, /* phi_1 and phi_2 != 0 */
  69                 phi_1_is_zero =  1, /* phi_1 = 0 */
  70                 phi_2_is_zero = -1  /* phi_2 = 0 */
  71             };
  72
  73             template <typename T>
  74             struct par_imw_p
  75             {
  76                 T    P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2;
  77                 T    phi_1, phi_2, lam_1;
  78                 detail::en<T> en;
  79                 mode_type mode;
  80             };
  81
  82             template <typename Params, typename T>
  83             inline int phi12(Params const& params,
  84                              par_imw_p<T> & proj_parm, T *del, T *sig)
  85             {
  86                 int err = 0;
  87
  88                 if (!pj_param_r<srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1, proj_parm.phi_1) ||
  89                     !pj_param_r<srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2, proj_parm.phi_2)) {
  90                     err = -41;
  91                 } else {
  92                     //proj_parm.phi_1 = pj_get_param_r(par.params, "lat_1"); // set above
  93                     //proj_parm.phi_2 = pj_get_param_r(par.params, "lat_2"); // set above
  94                     *del = 0.5 * (proj_parm.phi_2 - proj_parm.phi_1);
  95                     *sig = 0.5 * (proj_parm.phi_2 + proj_parm.phi_1);
  96                     err = (fabs(*del) < epsilon || fabs(*sig) < epsilon) ? -42 : 0;
  97                 }
  98                 return err;
  99             }
 100             template <typename Parameters, typename T>
 101             inline point_xy<T> loc_for(T const& lp_lam, T const& lp_phi,
 102                                        Parameters const& par,
 103                                        par_imw_p<T> const& proj_parm,
 104                                        T *yc)
 105             {
 106                 point_xy<T> xy;
 107
 108                 if (lp_phi == 0.0) {
 109                     xy.x = lp_lam;
 110                     xy.y = 0.;
 111                 } else {
 112                     T xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
 113
 114                     sp = sin(lp_phi);
 115                     m = pj_mlfn(lp_phi, sp, cos(lp_phi), proj_parm.en);
 116                     xa = proj_parm.Pp + proj_parm.Qp * m;
 117                     ya = proj_parm.P + proj_parm.Q * m;
 118                     R = 1. / (tan(lp_phi) * sqrt(1. - par.es * sp * sp));
 119                     C = sqrt(R * R - xa * xa);
 120                     if (lp_phi < 0.) C = - C;
 121                     C += ya - R;
 122                     if (proj_parm.mode == phi_2_is_zero) {
 123                         xb = lp_lam;
 124                         yb = proj_parm.C2;
 125                     } else {
 126                         t = lp_lam * proj_parm.sphi_2;
 127                         xb = proj_parm.R_2 * sin(t);
 128                         yb = proj_parm.C2 + proj_parm.R_2 * (1. - cos(t));
 129                     }
 130                     if (proj_parm.mode == phi_1_is_zero) {
 131                         xc = lp_lam;
 132                         *yc = 0.;
 133                     } else {
 134                         t = lp_lam * proj_parm.sphi_1;
 135                         xc = proj_parm.R_1 * sin(t);
 136                         *yc = proj_parm.R_1 * (1. - cos(t));
 137                     }
 138                     D = (xb - xc)/(yb - *yc);
 139                     B = xc + D * (C + R - *yc);
 140                     xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
 141                     if (lp_phi > 0)
 142                         xy.x = - xy.x;
 143                     xy.x = (B + xy.x) / (1. + D * D);
 144                     xy.y = sqrt(R * R - xy.x * xy.x);
 145                     if (lp_phi > 0)
 146                         xy.y = - xy.y;
 147                     xy.y += C + R;
 148                 }
 149                 return (xy);
 150             }
 151             template <typename Parameters, typename T>
 152             inline void xy(Parameters const& par, par_imw_p<T> const& proj_parm,
 153                            T const& phi,
 154                            T *x, T *y, T *sp, T *R)
 155             {
 156                 T F;
 157
 158                 *sp = sin(phi);
 159                 *R = 1./(tan(phi) * sqrt(1. - par.es * *sp * *sp ));
 160                 F = proj_parm.lam_1 * *sp;
 161                 *y = *R * (1 - cos(F));
 162                 *x = *R * sin(F);
 163             }
 164
 165             template <typename T, typename Parameters>
 166             struct base_imw_p_ellipsoid
 167             {
 168                 par_imw_p<T> m_proj_parm;
 169
 170                 // FORWARD(e_forward)  ellipsoid
 171                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 172                 inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 173                 {
 174                     T yc = 0;
 175                     point_xy<T> xy = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
 176                     xy_x = xy.x; xy_y = xy.y;
 177                 }
 178
 179                 // INVERSE(e_inverse)  ellipsoid
 180                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 181                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
 182                 {
 183                     point_xy<T> t;
 184                     T yc = 0.0;
 185                     int i = 0;
 186                     const int n_max_iter = 1000; /* Arbitrarily choosen number... */
 187
 188                     lp_lat = this->m_proj_parm.phi_2;
 189                     lp_lon = xy_x / cos(lp_lat);
 190                     do {
 191                         t = loc_for(lp_lon, lp_lat, par, m_proj_parm, &yc);
 192                         lp_lat = ((lp_lat - this->m_proj_parm.phi_1) * (xy_y - yc) / (t.y - yc)) + this->m_proj_parm.phi_1;
 193                         lp_lon = lp_lon * xy_x / t.x;
 194                         i++;
 195                     } while (i < n_max_iter &&
 196                              (fabs(t.x - xy_x) > tolerance || fabs(t.y - xy_y) > tolerance));
 197
 198                     if( i == n_max_iter )
 199                     {
 200                         lp_lon = lp_lat = HUGE_VAL;
 201                     }
 202                 }
 203
 204                 static inline std::string get_name()
 205                 {
 206                     return "imw_p_ellipsoid";
 207                 }
 208
 209             };
 210
 211             // International Map of the World Polyconic
 212             template <typename Params, typename Parameters, typename T>
 213             inline void setup_imw_p(Params const& params, Parameters const& par, par_imw_p<T>& proj_parm)
 214             {
 215                 T del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
 216                 int err;
 217
 218                 proj_parm.en = pj_enfn<T>(par.es);
 219                 if( (err = phi12(params, proj_parm, &del, &sig)) != 0)
 220                     BOOST_THROW_EXCEPTION( projection_exception(err) );
 221                 if (proj_parm.phi_2 < proj_parm.phi_1) { /* make sure proj_parm.phi_1 most southerly */
 222                     del = proj_parm.phi_1;
 223                     proj_parm.phi_1 = proj_parm.phi_2;
 224                     proj_parm.phi_2 = del;
 225                 }
 226                 if (pj_param_r<srs::spar::lon_1>(params, "lon_1", srs::dpar::lon_1, proj_parm.lam_1)) {
 227                     /* empty */
 228                 } else { /* use predefined based upon latitude */
 229                     sig = fabs(sig * geometry::math::r2d<T>());
 230                     if (sig <= 60)      sig = 2.;
 231                     else if (sig <= 76) sig = 4.;
 232                     else                sig = 8.;
 233                     proj_parm.lam_1 = sig * geometry::math::d2r<T>();
 234                 }
 235                 proj_parm.mode = none_is_zero;
 236                 if (proj_parm.phi_1 != 0.0)
 237                     xy(par, proj_parm, proj_parm.phi_1, &x1, &y1, &proj_parm.sphi_1, &proj_parm.R_1);
 238                 else {
 239                     proj_parm.mode = phi_1_is_zero;
 240                     y1 = 0.;
 241                     x1 = proj_parm.lam_1;
 242                 }
 243                 if (proj_parm.phi_2 != 0.0)
 244                     xy(par, proj_parm, proj_parm.phi_2, &x2, &T2, &proj_parm.sphi_2, &proj_parm.R_2);
 245                 else {
 246                     proj_parm.mode = phi_2_is_zero;
 247                     T2 = 0.;
 248                     x2 = proj_parm.lam_1;
 249                 }
 250                 m1 = pj_mlfn(proj_parm.phi_1, proj_parm.sphi_1, cos(proj_parm.phi_1), proj_parm.en);
 251                 m2 = pj_mlfn(proj_parm.phi_2, proj_parm.sphi_2, cos(proj_parm.phi_2), proj_parm.en);
 252                 t = m2 - m1;
 253                 s = x2 - x1;
 254                 y2 = sqrt(t * t - s * s) + y1;
 255                 proj_parm.C2 = y2 - T2;
 256                 t = 1. / t;
 257                 proj_parm.P = (m2 * y1 - m1 * y2) * t;
 258                 proj_parm.Q = (y2 - y1) * t;
 259                 proj_parm.Pp = (m2 * x1 - m1 * x2) * t;
 260                 proj_parm.Qp = (x2 - x1) * t;
 261             }
 262
 263     }} // namespace detail::imw_p
 264     #endif // doxygen
 265
 266     /*!
 267         \brief International Map of the World Polyconic projection
 268         \ingroup projections
 269         \tparam Geographic latlong point type
 270         \tparam Cartesian xy point type
 271         \tparam Parameters parameter type
 272         \par Projection characteristics
 273          - Mod. Polyconic
 274          - Ellipsoid
 275         \par Projection parameters
 276          - lat_1: Latitude of first standard parallel
 277          - lat_2: Latitude of second standard parallel
 278          - lon_1 (degrees)
 279         \par Example
 280         \image html ex_imw_p.gif
 281     */
 282     template <typename T, typename Parameters>
 283     struct imw_p_ellipsoid : public detail::imw_p::base_imw_p_ellipsoid<T, Parameters>
 284     {
 285         template <typename Params>
 286         inline imw_p_ellipsoid(Params const& params, Parameters const& par)
 287         {
 288             detail::imw_p::setup_imw_p(params, par, this->m_proj_parm);
 289         }
 290     };
 291
 292     #ifndef DOXYGEN_NO_DETAIL
 293     namespace detail
 294     {
 295
 296         // Static projection
 297         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_imw_p, imw_p_ellipsoid)
 298
 299         // Factory entry(s)
 300         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(imw_p_entry, imw_p_ellipsoid)
 301
 302         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(imw_p_init)
 303         {
 304             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(imw_p, imw_p_entry)
 305         }
 306
 307     } // namespace detail
 308     #endif // doxygen
 309
 310 } // namespace projections
 311
 312 }} // namespace boost::geometry
 313
 314 #endif // BOOST_GEOMETRY_PROJECTIONS_IMW_P_HPP
 315