projections/impl/pj_init.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002 // This file is manually converted from PROJ4
0003
0004 // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
0005 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0006
0007 // This file was modified by Oracle on 2017-2020.
0008 // Modifications copyright (c) 2017-2020, Oracle and/or its affiliates.
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010
0011 // Use, modification and distribution is subject to the Boost Software License,
0012 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0013 // http://www.boost.org/LICENSE_1_0.txt)
0014
0015 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0016 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0017 // PROJ4 is maintained by Frank Warmerdam
0018 // PROJ4 is converted to Geometry Library by Barend Gehrels (Geodan, Amsterdam)
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_IMPL_PJ_INIT_HPP
0041 #define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP
0042
0043 #include <cstdlib>
0044 #include <string>
0045 #include <type_traits>
0046 #include <vector>
0047
0048 #include <boost/geometry/core/static_assert.hpp>
0049
0050 #include <boost/geometry/srs/projections/dpar.hpp>
0051 #include <boost/geometry/srs/projections/impl/dms_parser.hpp>
0052 #include <boost/geometry/srs/projections/impl/pj_datum_set.hpp>
0053 #include <boost/geometry/srs/projections/impl/pj_datums.hpp>
0054 #include <boost/geometry/srs/projections/impl/pj_ell_set.hpp>
0055 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0056 #include <boost/geometry/srs/projections/impl/pj_units.hpp>
0057 #include <boost/geometry/srs/projections/impl/projects.hpp>
0058 #include <boost/geometry/srs/projections/proj4.hpp>
0059 #include <boost/geometry/srs/projections/spar.hpp>
0060
0061 #include <boost/geometry/util/condition.hpp>
0062 #include <boost/geometry/util/math.hpp>
0063
0064
0065 namespace boost { namespace geometry { namespace projections
0066 {
0067
0068
0069 namespace detail
0070 {
0071
0072 /************************************************************************/
0073 /*                           pj_init_proj()                             */
0074 /************************************************************************/
0075
0076 template <typename T>
0077 inline void pj_init_proj(srs::detail::proj4_parameters const& params,
0078                          parameters<T> & par)
0079 {
0080     par.id = pj_get_param_s(params, "proj");
0081 }
0082
0083 template <typename T>
0084 inline void pj_init_proj(srs::dpar::parameters<T> const& params,
0085                          parameters<T> & par)
0086 {
0087     typename srs::dpar::parameters<T>::const_iterator it = pj_param_find(params, srs::dpar::proj);
0088     if (it != params.end())
0089     {
0090         par.id = static_cast<srs::dpar::value_proj>(it->template get_value<int>());
0091     }
0092 }
0093
0094 template <typename T, typename ...Ps>
0095 inline void pj_init_proj(srs::spar::parameters<Ps...> const& ,
0096                          parameters<T> & par)
0097 {
0098     typedef srs::spar::parameters<Ps...> params_type;
0099     typedef typename geometry::tuples::find_if
0100         <
0101             params_type,
0102             srs::spar::detail::is_param_tr<srs::spar::detail::proj_traits>::pred
0103         >::type proj_type;
0104
0105     static const bool is_found = geometry::tuples::is_found<proj_type>::value;
0106
0107     BOOST_GEOMETRY_STATIC_ASSERT((is_found), "Projection not named.", params_type);
0108
0109     par.id = srs::spar::detail::proj_traits<proj_type>::id;
0110 }
0111
0112 /************************************************************************/
0113 /*                           pj_init_units()                            */
0114 /************************************************************************/
0115
0116 template <typename T, bool Vertical>
0117 inline void pj_init_units(srs::detail::proj4_parameters const& params,
0118                           T & to_meter,
0119                           T & fr_meter,
0120                           T const& default_to_meter,
0121                           T const& default_fr_meter)
0122 {
0123     std::string s;
0124     std::string units = pj_get_param_s(params, Vertical ? "vunits" : "units");
0125     if (! units.empty())
0126     {
0127         static const int n = sizeof(pj_units) / sizeof(pj_units[0]);
0128         int index = -1;
0129         for (int i = 0; i < n && index == -1; i++)
0130         {
0131             if(pj_units[i].id == units)
0132             {
0133                 index = i;
0134             }
0135         }
0136
0137         if (index == -1) {
0138             BOOST_THROW_EXCEPTION( projection_exception(error_unknow_unit_id) );
0139         }
0140         s = pj_units[index].to_meter;
0141     }
0142
0143     if (s.empty())
0144     {
0145         s = pj_get_param_s(params, Vertical ? "vto_meter" : "to_meter");
0146     }
0147
0148     // TODO: numerator and denominator could be taken from pj_units
0149     if (! s.empty())
0150     {
0151         std::size_t const pos = s.find('/');
0152         if (pos == std::string::npos)
0153         {
0154             to_meter = geometry::str_cast<T>(s);
0155         }
0156         else
0157         {
0158             T const numerator = geometry::str_cast<T>(s.substr(0, pos));
0159             T const denominator = geometry::str_cast<T>(s.substr(pos + 1));
0160             if (numerator == 0.0 || denominator == 0.0)
0161             {
0162                 BOOST_THROW_EXCEPTION( projection_exception(error_unit_factor_less_than_0) );
0163             }
0164             to_meter = numerator / denominator;
0165         }
0166         if (to_meter == 0.0)
0167         {
0168             BOOST_THROW_EXCEPTION( projection_exception(error_unit_factor_less_than_0) );
0169         }
0170         fr_meter = 1. / to_meter;
0171     }
0172     else
0173     {
0174         to_meter = default_to_meter;
0175         fr_meter = default_fr_meter;
0176     }
0177 }
0178
0179 template <typename T, bool Vertical>
0180 inline void pj_init_units(srs::dpar::parameters<T> const& params,
0181                           T & to_meter,
0182                           T & fr_meter,
0183                           T const& default_to_meter,
0184                           T const& default_fr_meter)
0185 {
0186     typename srs::dpar::parameters<T>::const_iterator
0187         it = pj_param_find(params, Vertical ? srs::dpar::vunits : srs::dpar::units);
0188     if (it != params.end())
0189     {
0190         static const int n = sizeof(pj_units) / sizeof(pj_units[0]);
0191         const int i = it->template get_value<int>();
0192         if (i >= 0 && i < n)
0193         {
0194             T const numerator = pj_units[i].numerator;
0195             T const denominator = pj_units[i].denominator;
0196             if (numerator == 0.0 || denominator == 0.0)
0197             {
0198                 BOOST_THROW_EXCEPTION( projection_exception(error_unit_factor_less_than_0) );
0199             }
0200             to_meter = numerator / denominator;
0201             fr_meter = 1. / to_meter;
0202         }
0203         else
0204         {
0205             BOOST_THROW_EXCEPTION( projection_exception(error_unknow_unit_id) );
0206         }
0207     }
0208     else
0209     {
0210         it = pj_param_find(params, Vertical ? srs::dpar::vto_meter : srs::dpar::to_meter);
0211         if (it != params.end())
0212         {
0213             to_meter = it->template get_value<T>();
0214             fr_meter = 1. / to_meter;
0215         }
0216         else
0217         {
0218             to_meter = default_to_meter;
0219             fr_meter = default_fr_meter;
0220         }
0221     }
0222 }
0223
0224 template
0225 <
0226     typename Params,
0227     bool Vertical,
0228     int UnitsI = geometry::tuples::find_index_if
0229         <
0230             Params,
0231             std::conditional_t
0232                 <
0233                     Vertical,
0234                     srs::spar::detail::is_param_t<srs::spar::vunits>,
0235                     srs::spar::detail::is_param_tr<srs::spar::detail::units_traits>
0236                 >::template pred
0237         >::value,
0238     int ToMeterI = geometry::tuples::find_index_if
0239         <
0240             Params,
0241             std::conditional_t
0242                 <
0243                     Vertical,
0244                     srs::spar::detail::is_param_t<srs::spar::vto_meter>,
0245                     srs::spar::detail::is_param_t<srs::spar::to_meter>
0246                 >::template pred
0247         >::value,
0248     int N = geometry::tuples::size<Params>::value
0249 >
0250 struct pj_init_units_static
0251     : pj_init_units_static<Params, Vertical, UnitsI, N, N>
0252 {};
0253
0254 template <typename Params, bool Vertical, int UnitsI, int N>
0255 struct pj_init_units_static<Params, Vertical, UnitsI, N, N>
0256 {
0257     static const int n = sizeof(pj_units) / sizeof(pj_units[0]);
0258     static const int i = srs::spar::detail::units_traits
0259                     <
0260                         typename geometry::tuples::element<UnitsI, Params>::type
0261                     >::id;
0262     static const bool is_valid = i >= 0 && i < n;
0263
0264     BOOST_GEOMETRY_STATIC_ASSERT((is_valid), "Unknown unit ID.", Params);
0265
0266     template <typename T>
0267     static void apply(Params const& ,
0268                       T & to_meter, T & fr_meter,
0269                       T const& , T const& )
0270     {
0271         T const numerator = pj_units[i].numerator;
0272         T const denominator = pj_units[i].denominator;
0273         if (numerator == 0.0 || denominator == 0.0)
0274         {
0275             BOOST_THROW_EXCEPTION( projection_exception(error_unit_factor_less_than_0) );
0276         }
0277         to_meter = numerator / denominator;
0278         fr_meter = 1. / to_meter;
0279     }
0280 };
0281
0282 template <typename Params, bool Vertical, int ToMeterI, int N>
0283 struct pj_init_units_static<Params, Vertical, N, ToMeterI, N>
0284 {
0285     template <typename T>
0286     static void apply(Params const& params,
0287                       T & to_meter, T & fr_meter,
0288                       T const& , T const& )
0289     {
0290         to_meter = geometry::tuples::get<ToMeterI>(params).value;
0291         fr_meter = 1. / to_meter;
0292     }
0293 };
0294
0295 template <typename Params, bool Vertical, int N>
0296 struct pj_init_units_static<Params, Vertical, N, N, N>
0297 {
0298     template <typename T>
0299     static void apply(Params const& ,
0300                       T & to_meter, T & fr_meter,
0301                       T const& default_to_meter, T const& default_fr_meter)
0302     {
0303         to_meter = default_to_meter;
0304         fr_meter = default_fr_meter;
0305     }
0306 };
0307
0308 template <typename T, bool Vertical, typename ...Ps>
0309 inline void pj_init_units(srs::spar::parameters<Ps...> const& params,
0310                           T & to_meter,
0311                           T & fr_meter,
0312                           T const& default_to_meter,
0313                           T const& default_fr_meter)
0314 {
0315     pj_init_units_static
0316         <
0317             srs::spar::parameters<Ps...>,
0318             Vertical
0319         >::apply(params, to_meter, fr_meter, default_to_meter, default_fr_meter);
0320 }
0321
0322 /************************************************************************/
0323 /*                           pj_init_pm()                               */
0324 /************************************************************************/
0325
0326 template <typename T>
0327 inline void pj_init_pm(srs::detail::proj4_parameters const& params, T& val)
0328 {
0329     std::string pm = pj_get_param_s(params, "pm");
0330     if (! pm.empty())
0331     {
0332         int n = sizeof(pj_prime_meridians) / sizeof(pj_prime_meridians[0]);
0333         for (int i = 0; i < n ; i++)
0334         {
0335             if(pj_prime_meridians[i].id == pm)
0336             {
0337                 val = pj_prime_meridians[i].deg * math::d2r<T>();
0338                 return;
0339             }
0340         }
0341
0342         // TODO: Is this try-catch needed?
0343         // In other cases the bad_str_cast exception is simply thrown
0344         BOOST_TRY
0345         {
0346             val = dms_parser<T, true>::apply(pm).angle();
0347             return;
0348         }
0349         BOOST_CATCH(geometry::bad_str_cast const&)
0350         {
0351             BOOST_THROW_EXCEPTION( projection_exception(error_unknown_prime_meridian) );
0352         }
0353         BOOST_CATCH_END
0354     }
0355
0356     val = 0.0;
0357 }
0358
0359 template <typename T>
0360 inline void pj_init_pm(srs::dpar::parameters<T> const& params, T& val)
0361 {
0362     typename srs::dpar::parameters<T>::const_iterator it = pj_param_find(params, srs::dpar::pm);
0363     if (it != params.end())
0364     {
0365         if (it->template is_value_set<int>())
0366         {
0367             int n = sizeof(pj_prime_meridians) / sizeof(pj_prime_meridians[0]);
0368             int i = it->template get_value<int>();
0369             if (i >= 0 && i < n)
0370             {
0371                 val = pj_prime_meridians[i].deg * math::d2r<T>();
0372                 return;
0373             }
0374             else
0375             {
0376                 BOOST_THROW_EXCEPTION( projection_exception(error_unknown_prime_meridian) );
0377             }
0378         }
0379         else if (it->template is_value_set<T>())
0380         {
0381             val = it->template get_value<T>() * math::d2r<T>();
0382             return;
0383         }
0384     }
0385
0386     val = 0.0;
0387 }
0388
0389 template
0390 <
0391     typename Params,
0392     int I = geometry::tuples::find_index_if
0393         <
0394             Params,
0395             srs::spar::detail::is_param_tr<srs::spar::detail::pm_traits>::pred
0396         >::value,
0397     int N = geometry::tuples::size<Params>::value
0398 >
0399 struct pj_init_pm_static
0400 {
0401     template <typename T>
0402     static void apply(Params const& params, T & val)
0403     {
0404         typedef typename geometry::tuples::element<I, Params>::type param_type;
0405
0406         val = srs::spar::detail::pm_traits<param_type>::value(geometry::tuples::get<I>(params));
0407     }
0408 };
0409 template <typename Params, int N>
0410 struct pj_init_pm_static<Params, N, N>
0411 {
0412     template <typename T>
0413     static void apply(Params const& , T & val)
0414     {
0415         val = 0;
0416     }
0417 };
0418
0419 template <typename T, typename ...Ps>
0420 inline void pj_init_pm(srs::spar::parameters<Ps...> const& params, T& val)
0421 {
0422     pj_init_pm_static
0423         <
0424             srs::spar::parameters<Ps...>
0425         >::apply(params, val);
0426 }
0427
0428 /************************************************************************/
0429 /*                        pj_init_axis()                                */
0430 /************************************************************************/
0431
0432 template <typename Params, typename T>
0433 inline void pj_init_axis(Params const& params, parameters<T> & projdef)
0434 {
0435     std::string axis = pj_get_param_s(params, "axis");
0436     if(! axis.empty())
0437     {
0438         for (std::size_t i = 0; i < axis.length(); ++i)
0439         {
0440             switch(axis[i])
0441             {
0442                 case 'w':
0443                     projdef.sign[i] = -1;
0444                     projdef.axis[i] = 0;
0445                     break;
0446                 case 'e':
0447                     projdef.sign[i] = 1;
0448                     projdef.axis[i] = 0;
0449                     break;
0450                 case 's':
0451                     projdef.sign[i] = -1;
0452                     projdef.axis[i] = 1;
0453                     break;
0454                 case 'n':
0455                     projdef.sign[i] = 1;
0456                     projdef.axis[i] = 1;
0457                     break;
0458                 case 'd':
0459                     projdef.sign[i] = -1;
0460                     projdef.axis[i] = 2;
0461                     break;
0462                 case 'u':
0463                     projdef.sign[i] = 1;
0464                     projdef.axis[i] = 2;
0465                     break;
0466                 default:
0467                     BOOST_THROW_EXCEPTION( projection_exception(error_axis) );
0468             }
0469         }
0470         // Currently not support elevation
0471         if (projdef.axis[0] + projdef.axis[1] != 1)
0472         {
0473             BOOST_THROW_EXCEPTION( projection_exception(error_axis) );
0474         }
0475     }
0476
0477 }
0478
0479 // TODO: implement axis support for other types of parameters
0480
0481 template <typename T>
0482 inline void pj_init_axis(srs::dpar::parameters<T> const& , parameters<T> & )
0483 {}
0484
0485 template <typename Params>
0486 struct pj_init_axis_static
0487 {
0488     template <typename T>
0489     static void apply(Params const& , parameters<T> & )
0490     {}
0491 };
0492
0493 template <typename T, typename ...Ps>
0494 inline void pj_init_axis(srs::spar::parameters<Ps...> const& params, parameters<T> & projdef)
0495 {
0496     pj_init_axis_static
0497         <
0498             srs::spar::parameters<Ps...>
0499         >::apply(params, projdef);
0500 }
0501
0502
0503 /************************************************************************/
0504 /*                              pj_init()                               */
0505 /*                                                                      */
0506 /*      Main entry point for initialing a PJ projections                */
0507 /*      definition.                                                     */
0508 /************************************************************************/
0509 template <typename T, typename Params>
0510 inline parameters<T> pj_init(Params const& params)
0511 {
0512     parameters<T> pin;
0513
0514     // find projection -> implemented in projection factory
0515     pj_init_proj(params, pin);
0516     // exception thrown in projection<>
0517     // TODO: consider throwing here both projection_unknown_id_exception and
0518     // projection_not_named_exception in order to throw before other exceptions
0519     //if (pin.name.empty())
0520     //{ BOOST_THROW_EXCEPTION( projection_not_named_exception() ); }
0521
0522     // NOTE: proj4 gets defaults from "proj_def.dat".
0523     // In Boost.Geometry this is emulated by manually setting them in
0524     // pj_ell_init and projections aea, lcc and lagrng
0525
0526     /* set datum parameters */
0527     pj_datum_init(params, pin);
0528
0529     /* set ellipsoid/sphere parameters */
0530     pj_ell_init(params, pin.a, pin.es);
0531
0532     pin.a_orig = pin.a;
0533     pin.es_orig = pin.es;
0534
0535     pin.e = sqrt(pin.es);
0536     pin.ra = 1. / pin.a;
0537     pin.one_es = 1. - pin.es;
0538     if (pin.one_es == 0.) {
0539         BOOST_THROW_EXCEPTION( projection_exception(error_eccentricity_is_one) );
0540     }
0541     pin.rone_es = 1./pin.one_es;
0542
0543     /* Now that we have ellipse information check for WGS84 datum */
0544     if( pin.datum_type == datum_3param
0545         && pin.datum_params[0] == 0.0
0546         && pin.datum_params[1] == 0.0
0547         && pin.datum_params[2] == 0.0
0548         && pin.a == 6378137.0
0549         && geometry::math::abs(pin.es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/
0550     {
0551         pin.datum_type = datum_wgs84;
0552     }
0553
0554     /* set pin.geoc coordinate system */
0555     pin.geoc = (pin.es && pj_get_param_b<srs::spar::geoc>(params, "geoc", srs::dpar::geoc));
0556
0557     /* over-ranging flag */
0558     pin.over = pj_get_param_b<srs::spar::over>(params, "over", srs::dpar::over);
0559
0560     /* longitude center for wrapping */
0561     pin.is_long_wrap_set = pj_param_r<srs::spar::lon_wrap>(params, "lon_wrap", srs::dpar::lon_wrap, pin.long_wrap_center);
0562
0563     /* central meridian */
0564     pin.lam0 = pj_get_param_r<T, srs::spar::lon_0>(params, "lon_0", srs::dpar::lon_0);
0565
0566     /* central latitude */
0567     pin.phi0 = pj_get_param_r<T, srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0);
0568
0569     /* false easting and northing */
0570     pin.x0 = pj_get_param_f<T, srs::spar::x_0>(params, "x_0", srs::dpar::x_0);
0571     pin.y0 = pj_get_param_f<T, srs::spar::y_0>(params, "y_0", srs::dpar::y_0);
0572
0573     /* general scaling factor */
0574     if (pj_param_f<srs::spar::k_0>(params, "k_0", srs::dpar::k_0, pin.k0)) {
0575         /* empty */
0576     } else if (pj_param_f<srs::spar::k>(params, "k", srs::dpar::k, pin.k0)) {
0577         /* empty */
0578     } else
0579         pin.k0 = 1.;
0580     if (pin.k0 <= 0.) {
0581         BOOST_THROW_EXCEPTION( projection_exception(error_k_less_than_zero) );
0582     }
0583
0584     /* set units */
0585     pj_init_units<T, false>(params, pin.to_meter, pin.fr_meter, 1., 1.);
0586     pj_init_units<T, true>(params, pin.vto_meter, pin.vfr_meter, pin.to_meter, pin.fr_meter);
0587
0588     /* prime meridian */
0589     pj_init_pm(params, pin.from_greenwich);
0590
0591     /* set axis orientation */
0592     pj_init_axis(params, pin);
0593
0594     return pin;
0595 }
0596
0597
0598 } // namespace detail
0599 }}} // namespace boost::geometry::projections
0600
0601 #endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP