projections/impl/pj_ell_set.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
0006 // This file was modified by Oracle on 2017-2020.
0007 // Modifications copyright (c) 2017-2020, Oracle and/or its affiliates.
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 Geometry Library by
0018 //     Barend Gehrels (Geodan, Amsterdam)
0019 //     Adam Wulkiewicz
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_IMPL_PJ_ELL_SET_HPP
0042 #define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP
0043
0044 #include <string>
0045 #include <type_traits>
0046 #include <vector>
0047
0048 #include <boost/geometry/core/static_assert.hpp>
0049 #include <boost/geometry/formulas/eccentricity_sqr.hpp>
0050 #include <boost/geometry/util/math.hpp>
0051
0052 #include <boost/geometry/srs/projections/dpar.hpp>
0053 #include <boost/geometry/srs/projections/impl/pj_datum_set.hpp>
0054 #include <boost/geometry/srs/projections/impl/pj_ellps.hpp>
0055 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0056 #include <boost/geometry/srs/projections/proj4.hpp>
0057 #include <boost/geometry/srs/projections/spar.hpp>
0058
0059
0060 namespace boost { namespace geometry { namespace projections {
0061
0062 namespace detail {
0063
0064 /* set ellipsoid parameters a and es */
0065 template <typename T>
0066 inline T SIXTH() { return .1666666666666666667; } /* 1/6 */
0067 template <typename T>
0068 inline T RA4() { return .04722222222222222222; } /* 17/360 */
0069 template <typename T>
0070 inline T RA6() { return .02215608465608465608; } /* 67/3024 */
0071 template <typename T>
0072 inline T RV4() { return .06944444444444444444; } /* 5/72 */
0073 template <typename T>
0074 inline T RV6() { return .04243827160493827160; } /* 55/1296 */
0075
0076 template <typename T>
0077 inline T pj_ell_b_to_es(T const& a, T const& b)
0078 {
0079     return 1. - (b * b) / (a * a);
0080 }
0081
0082 /************************************************************************/
0083 /*                          pj_ell_init_ellps()                         */
0084 /************************************************************************/
0085
0086 // Originally a part of pj_ell_set()
0087 template <typename T>
0088 inline bool pj_ell_init_ellps(srs::detail::proj4_parameters const& params, T &a, T &b)
0089 {
0090     /* check if ellps present and temporarily append its values to pl */
0091     std::string name = pj_get_param_s(params, "ellps");
0092     if (! name.empty())
0093     {
0094         const pj_ellps_type<T>* pj_ellps = pj_get_ellps<T>().first;
0095         const int n = pj_get_ellps<T>().second;
0096         int index = -1;
0097         for (int i = 0; i < n && index == -1; i++)
0098         {
0099             if(pj_ellps[i].id == name)
0100             {
0101                 index = i;
0102             }
0103         }
0104
0105         if (index == -1) {
0106             BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) );
0107         }
0108
0109         pj_ellps_type<T> const& pj_ellp = pj_ellps[index];
0110         a = pj_ellp.a;
0111         b = pj_ellp.b;
0112
0113         return true;
0114     }
0115
0116     return false;
0117 }
0118
0119 template <typename T>
0120 inline bool pj_ell_init_ellps(srs::dpar::parameters<T> const& params, T &a, T &b)
0121 {
0122     /* check if ellps present and temporarily append its values to pl */
0123     typename srs::dpar::parameters<T>::const_iterator
0124         it = pj_param_find(params, srs::dpar::ellps);
0125     if (it != params.end())
0126     {
0127         if (it->template is_value_set<int>())
0128         {
0129             const pj_ellps_type<T>* pj_ellps = pj_get_ellps<T>().first;
0130             const int n = pj_get_ellps<T>().second;
0131             int i = it->template get_value<int>();
0132
0133             if (i < 0 || i >= n) {
0134                 BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) );
0135             }
0136
0137             pj_ellps_type<T> const& pj_ellp = pj_ellps[i];
0138             a = pj_ellp.a;
0139             b = pj_ellp.b;
0140         }
0141         else if (it->template is_value_set<T>())
0142         {
0143             a = it->template get_value<T>();
0144             b = a;
0145         }
0146         else if (it->template is_value_set<srs::spheroid<T> >())
0147         {
0148             srs::spheroid<T> const& s = it->template get_value<srs::spheroid<T> >();
0149             a = geometry::get_radius<0>(s);
0150             b = geometry::get_radius<2>(s);
0151         }
0152         else
0153         {
0154             BOOST_THROW_EXCEPTION( projection_exception(error_unknown_ellp_param) );
0155         }
0156
0157         return true;
0158     }
0159
0160     return false;
0161 }
0162
0163 template
0164 <
0165     typename Params,
0166     int I = geometry::tuples::find_index_if
0167         <
0168             Params,
0169             srs::spar::detail::is_param_tr<srs::spar::detail::ellps_traits>::pred
0170         >::value,
0171     int N = geometry::tuples::size<Params>::value
0172 >
0173 struct pj_ell_init_ellps_static
0174 {
0175     template <typename T>
0176     static bool apply(Params const& params, T &a, T &b)
0177     {
0178         typedef typename geometry::tuples::element<I, Params>::type param_type;
0179         typedef srs::spar::detail::ellps_traits<param_type> traits_type;
0180         typedef typename traits_type::template model_type<T>::type model_type;
0181
0182         param_type const& param = geometry::tuples::get<I>(params);
0183         model_type const& model = traits_type::template model<T>(param);
0184
0185         a = geometry::get_radius<0>(model);
0186         b = geometry::get_radius<2>(model);
0187
0188         return true;
0189     }
0190 };
0191 template <typename Params, int N>
0192 struct pj_ell_init_ellps_static<Params, N, N>
0193 {
0194     template <typename T>
0195     static bool apply(Params const& , T & , T & )
0196     {
0197         return false;
0198     }
0199 };
0200
0201 template <typename T, typename ...Ps>
0202 inline bool pj_ell_init_ellps(srs::spar::parameters<Ps...> const& params,
0203                               T &a, T &b)
0204 {
0205     return pj_ell_init_ellps_static
0206         <
0207             srs::spar::parameters<Ps...>
0208         >::apply(params, a, b);
0209 }
0210
0211 /************************************************************************/
0212 /*                             pj_ell_init()                            */
0213 /************************************************************************/
0214
0215 /* initialize geographic shape parameters */
0216 // This function works differently than the original pj_ell_set().
0217 // It doesn't push parameters defined in ellps into params list.
0218 // Instead it tries to use size (a, R) and shape (es, e, rf, f, b) parameters
0219 // and then if needed falls back to ellps, then to datum and then to the default WGS84
0220 template <typename Params, typename T>
0221 inline void pj_ell_init(Params const& params, T &a, T &es)
0222 {
0223     /* check for varying forms of ellipsoid input */
0224     a = es = 0.;
0225
0226     /* R takes precedence */
0227     if (pj_param_f<srs::spar::r>(params, "R", srs::dpar::r, a)) {
0228         /* empty */
0229     } else { /* probable elliptical figure */
0230
0231         // Set ellipsoid's size parameter
0232         a = pj_get_param_f<T, srs::spar::a>(params, "a", srs::dpar::a);
0233         bool is_a_set = a != 0.0;
0234
0235         // Set ellipsoid's shape parameter
0236         T b = 0.0;
0237         bool is_ell_set = false;
0238         if (pj_param_f<srs::spar::es>(params, "es", srs::dpar::es, es)) {/* eccentricity squared */
0239             /* empty */
0240             is_ell_set = true;
0241         } else if (pj_param_f<srs::spar::e>(params, "e", srs::dpar::e, es)) { /* eccentricity */
0242             es = es * es;
0243             is_ell_set = true;
0244         } else if (pj_param_f<srs::spar::rf>(params, "rf", srs::dpar::rf, es)) { /* recip flattening */
0245             if (es == 0.0) {
0246                 BOOST_THROW_EXCEPTION( projection_exception(error_rev_flattening_is_zero) );
0247             }
0248             es = 1./ es;
0249             es = es * (2. - es);
0250             is_ell_set = true;
0251         } else if (pj_param_f<srs::spar::f>(params, "f", srs::dpar::f, es)) { /* flattening */
0252             es = es * (2. - es);
0253             is_ell_set = true;
0254         } else if (pj_param_f<srs::spar::b>(params, "b", srs::dpar::b, b)) { /* minor axis */
0255             es = pj_ell_b_to_es(a, b);
0256             is_ell_set = true;
0257         } /* else es == 0. and sphere of radius a */
0258
0259         // NOTE: Below when ellps is used to initialize a and es
0260         // b is not set because it only has sense together with a
0261         // but a could have been set separately before, e.g. consider passing:
0262         // a=1 ellps=airy (a=6377563.396 b=6356256.910)
0263         // after setting size parameter a and shape parameter from ellps
0264         // b has to be recalculated
0265
0266         // If ellipsoid's parameters are not set directly
0267         //   use ellps parameter
0268         if (! is_a_set || ! is_ell_set) {
0269             T ellps_a = 0, ellps_b = 0;
0270             if (pj_ell_init_ellps(params, ellps_a, ellps_b)) {
0271                 if (! is_a_set) {
0272                     a = ellps_a;
0273                     is_a_set = true;
0274                 }
0275                 if (! is_ell_set) {
0276                     es = pj_ell_b_to_es(ellps_a, ellps_b);
0277                     is_ell_set = true;
0278                 }
0279             }
0280         }
0281
0282         // If ellipsoid's parameters are not set
0283         //   use ellps defined by datum parameter
0284         if (! is_a_set || ! is_ell_set)
0285         {
0286             const pj_datums_type<T>* datum = pj_datum_find_datum<T>(params);
0287             if (datum != NULL)
0288             {
0289                 pj_ellps_type<T> const& pj_ellp = pj_get_ellps<T>().first[datum->ellps];
0290                 if (! is_a_set) {
0291                     a = pj_ellp.a;
0292                     is_a_set = true;
0293                 }
0294                 if (! is_ell_set) {
0295                     es = pj_ell_b_to_es(pj_ellp.a, pj_ellp.b);
0296                     is_ell_set = true;
0297                 }
0298             }
0299         }
0300
0301         // If ellipsoid's parameters are still not set
0302         //   use default WGS84
0303         if ((! is_a_set || ! is_ell_set)
0304          && ! pj_get_param_b<srs::spar::no_defs>(params, "no_defs", srs::dpar::no_defs))
0305         {
0306             pj_ellps_type<T> const& pj_ellp = pj_get_ellps<T>().first[srs::dpar::ellps_wgs84];
0307             if (! is_a_set) {
0308                 a = pj_ellp.a;
0309                 is_a_set = true;
0310             }
0311             if (! is_ell_set) {
0312                 es = pj_ell_b_to_es(pj_ellp.a, pj_ellp.b);
0313                 is_ell_set = true;
0314             }
0315         }
0316
0317         if (b == 0.0)
0318             b = a * sqrt(1. - es);
0319
0320         /* following options turn ellipsoid into equivalent sphere */
0321         if (pj_get_param_b<srs::spar::r_au>(params, "R_A", srs::dpar::r_au)) { /* sphere--area of ellipsoid */
0322             a *= 1. - es * (SIXTH<T>() + es * (RA4<T>() + es * RA6<T>()));
0323             es = 0.;
0324         } else if (pj_get_param_b<srs::spar::r_v>(params, "R_V", srs::dpar::r_v)) { /* sphere--vol. of ellipsoid */
0325             a *= 1. - es * (SIXTH<T>() + es * (RV4<T>() + es * RV6<T>()));
0326             es = 0.;
0327         } else if (pj_get_param_b<srs::spar::r_a>(params, "R_a", srs::dpar::r_a)) { /* sphere--arithmetic mean */
0328             a = .5 * (a + b);
0329             es = 0.;
0330         } else if (pj_get_param_b<srs::spar::r_g>(params, "R_g", srs::dpar::r_g)) { /* sphere--geometric mean */
0331             a = sqrt(a * b);
0332             es = 0.;
0333         } else if (pj_get_param_b<srs::spar::r_h>(params, "R_h", srs::dpar::r_h)) { /* sphere--harmonic mean */
0334             a = 2. * a * b / (a + b);
0335             es = 0.;
0336         } else {
0337             T tmp;
0338             bool i = pj_param_r<srs::spar::r_lat_a>(params, "R_lat_a", srs::dpar::r_lat_a, tmp);
0339             if (i || /* sphere--arith. */
0340                 pj_param_r<srs::spar::r_lat_g>(params, "R_lat_g", srs::dpar::r_lat_g, tmp)) { /* or geom. mean at latitude */
0341
0342                 tmp = sin(tmp);
0343                 if (geometry::math::abs(tmp) > geometry::math::half_pi<T>()) {
0344                     BOOST_THROW_EXCEPTION( projection_exception(error_ref_rad_larger_than_90) );
0345                 }
0346                 tmp = 1. - es * tmp * tmp;
0347                 a *= i ? .5 * (1. - es + tmp) / ( tmp * sqrt(tmp)) :
0348                     sqrt(1. - es) / tmp;
0349                 es = 0.;
0350             }
0351         }
0352     }
0353
0354     /* some remaining checks */
0355     if (es < 0.) {
0356         BOOST_THROW_EXCEPTION( projection_exception(error_es_less_than_zero) );
0357     }
0358     if (a <= 0.) {
0359         BOOST_THROW_EXCEPTION( projection_exception(error_major_axis_not_given) );
0360     }
0361 }
0362
0363 template <typename Params>
0364 struct static_srs_tag_check_nonexpanded
0365 {
0366     typedef std::conditional_t
0367         <
0368             geometry::tuples::exists_if
0369                 <
0370                     Params, srs::spar::detail::is_param_t<srs::spar::r>::pred
0371                 >::value
0372          || geometry::tuples::exists_if
0373                 <
0374                     Params, srs::spar::detail::is_param<srs::spar::r_au>::pred
0375                 >::value
0376          || geometry::tuples::exists_if
0377                 <
0378                     Params, srs::spar::detail::is_param<srs::spar::r_v>::pred
0379                 >::value
0380          || geometry::tuples::exists_if
0381                 <
0382                     Params, srs::spar::detail::is_param<srs::spar::r_a>::pred
0383                 >::value
0384          || geometry::tuples::exists_if
0385                 <
0386                     Params, srs::spar::detail::is_param<srs::spar::r_g>::pred
0387                 >::value
0388          || geometry::tuples::exists_if
0389                 <
0390                     Params, srs::spar::detail::is_param<srs::spar::r_h>::pred
0391                 >::value
0392          || geometry::tuples::exists_if
0393                 <
0394                     Params, srs::spar::detail::is_param_t<srs::spar::r_lat_a>::pred
0395                 >::value
0396          || geometry::tuples::exists_if
0397                 <
0398                     Params, srs::spar::detail::is_param_t<srs::spar::r_lat_g>::pred
0399                 >::value,
0400             srs_sphere_tag,
0401             // NOTE: The assumption here is that if the user defines either one of:
0402             // b, es, e, f, rf parameters then he wants to define spheroid, not sphere
0403             std::conditional_t
0404                 <
0405                     geometry::tuples::exists_if
0406                         <
0407                             Params, srs::spar::detail::is_param_t<srs::spar::b>::pred
0408                         >::value
0409                  || geometry::tuples::exists_if
0410                         <
0411                             Params, srs::spar::detail::is_param_t<srs::spar::es>::pred
0412                         >::value
0413                  || geometry::tuples::exists_if
0414                         <
0415                             Params, srs::spar::detail::is_param_t<srs::spar::e>::pred
0416                         >::value
0417                  || geometry::tuples::exists_if
0418                         <
0419                             Params, srs::spar::detail::is_param_t<srs::spar::rf>::pred
0420                         >::value
0421                  || geometry::tuples::exists_if
0422                         <
0423                             Params, srs::spar::detail::is_param_t<srs::spar::f>::pred
0424                         >::value,
0425                     srs_spheroid_tag,
0426                     void
0427                 >
0428         > type;
0429 };
0430
0431 template <typename Params>
0432 struct static_srs_tag_check_ellps
0433 {
0434     typedef typename geometry::tag
0435         <
0436             typename srs::spar::detail::ellps_traits
0437                 <
0438                     typename geometry::tuples::find_if
0439                         <
0440                             Params,
0441                             srs::spar::detail::is_param_tr<srs::spar::detail::ellps_traits>::pred
0442                         >::type
0443                 >::template model_type<double>::type // dummy type
0444         >::type type;
0445 };
0446
0447 template <typename Params>
0448 struct static_srs_tag_check_datum
0449 {
0450     typedef typename geometry::tag
0451         <
0452             typename srs::spar::detail::ellps_traits
0453                 <
0454                     typename srs::spar::detail::datum_traits
0455                         <
0456                             typename geometry::tuples::find_if
0457                                 <
0458                                     Params,
0459                                     srs::spar::detail::is_param_tr<srs::spar::detail::datum_traits>::pred
0460                                 >::type
0461                         >::ellps_type
0462                 >::template model_type<double>::type // dummy type
0463         >::type type;
0464 };
0465
0466 template
0467 <
0468     typename Params,
0469     typename NonExpandedTag = typename static_srs_tag_check_nonexpanded
0470                                 <
0471                                     Params
0472                                 >::type,
0473     typename EllpsTag = typename static_srs_tag_check_ellps
0474                             <
0475                                 Params
0476                             >::type,
0477     typename DatumTag = typename static_srs_tag_check_datum
0478                             <
0479                                 Params
0480                             >::type
0481 >
0482 struct static_srs_tag
0483 {
0484     // User passed one of the non-ellps, non-datum parameters
0485     typedef NonExpandedTag type;
0486 };
0487
0488 template <typename Params, typename EllpsTag, typename DatumTag>
0489 struct static_srs_tag<Params, void, EllpsTag, DatumTag>
0490 {
0491     // User didn't pass neither one of the non-ellps, non-datum parameters
0492     // but passed ellps
0493     typedef EllpsTag type;
0494 };
0495
0496 template <typename Params, typename DatumTag>
0497 struct static_srs_tag<Params, void, void, DatumTag>
0498 {
0499     // User didn't pass neither one of the non-ellps, non-datum parameters
0500     // nor ellps parameter but passed datum parameter
0501     typedef DatumTag type;
0502 };
0503
0504 template <typename Params>
0505 struct static_srs_tag<Params, void, void, void>
0506 {
0507     // User didn't pass any parameter defining model
0508     // so use default or generate error
0509     typedef std::conditional_t
0510         <
0511             geometry::tuples::exists_if
0512                 <
0513                     Params, srs::spar::detail::is_param<srs::spar::no_defs>::pred
0514                 >::value,
0515             void,
0516             srs_spheroid_tag // WGS84
0517         > type;
0518
0519     static const bool is_found = ! std::is_void<type>::value;
0520     BOOST_GEOMETRY_STATIC_ASSERT((is_found),
0521         "Expected ellipsoid or sphere definition.",
0522         Params);
0523 };
0524
0525
0526 template <typename T>
0527 inline void pj_calc_ellipsoid_params(parameters<T> & p, T const& a, T const& es) {
0528 /****************************************************************************************
0529     Calculate a large number of ancillary ellipsoidal parameters, in addition to
0530     the two traditional PROJ defining parameters: Semimajor axis, a, and the
0531     eccentricity squared, es.
0532
0533     Most of these parameters are fairly cheap to compute in comparison to the overall
0534     effort involved in initializing a PJ object. They may, however, take a substantial
0535     part of the time taken in computing an individual point transformation.
0536
0537     So by providing them up front, we can amortize the (already modest) cost over all
0538     transformations carried out over the entire lifetime of a PJ object, rather than
0539     incur that cost for every single transformation.
0540
0541     Most of the parameter calculations here are based on the "angular eccentricity",
0542     i.e. the angle, measured from the semiminor axis, of a line going from the north
0543     pole to one of the foci of the ellipsoid - or in other words: The arc sine of the
0544     eccentricity.
0545
0546     The formulae used are mostly taken from:
0547
0548     Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991).
0549     Columbus, Ohio:  Dept. of Geodetic Science
0550     and Surveying, Ohio State University.
0551
0552 ****************************************************************************************/
0553
0554     p.a = a;
0555     p.es = es;
0556
0557     /* Compute some ancillary ellipsoidal parameters */
0558     if (p.e==0)
0559         p.e = sqrt(p.es);  /* eccentricity */
0560     //p.alpha = asin (p.e);  /* angular eccentricity */
0561
0562     /* second eccentricity */
0563     //p.e2  = tan (p.alpha);
0564     //p.e2s = p.e2 * p.e2;
0565
0566     /* third eccentricity */
0567     //p.e3    = (0!=p.alpha)? sin (p.alpha) / sqrt(2 - sin (p.alpha)*sin (p.alpha)): 0;
0568     //p.e3s = p.e3 * p.e3;
0569
0570     /* flattening */
0571     //if (0==p.f)
0572     //    p.f  = 1 - cos (p.alpha);   /* = 1 - sqrt (1 - PIN->es); */
0573     //p.rf = p.f != 0.0 ? 1.0/p.f: HUGE_VAL;
0574
0575     /* second flattening */
0576     //p.f2  = (cos(p.alpha)!=0)? 1/cos (p.alpha) - 1: 0;
0577     //p.rf2 = p.f2 != 0.0 ? 1/p.f2: HUGE_VAL;
0578
0579     /* third flattening */
0580     //p.n  = pow (tan (p.alpha/2), 2);
0581     //p.rn = p.n != 0.0 ? 1/p.n: HUGE_VAL;
0582
0583     /* ...and a few more */
0584     //if (0==p.b)
0585     //    p.b  = (1 - p.f)*p.a;
0586     //p.rb = 1. / p.b;
0587     p.ra = 1. / p.a;
0588
0589     p.one_es = 1. - p.es;
0590     if (p.one_es == 0.) {
0591         BOOST_THROW_EXCEPTION( projection_exception(error_eccentricity_is_one) );
0592     }
0593
0594     p.rone_es = 1./p.one_es;
0595 }
0596
0597
0598 } // namespace detail
0599 }}} // namespace boost::geometry::projections
0600
0601 #endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_ELL_SET_HPP