ceph/src/boost/boost/geometry/srs/projections/impl/pj_init.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2 // This file is manually converted from PROJ4
   3
   4 // Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands.
   5
   6 // This file was modified by Oracle on 2017, 2018.
   7 // Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
   8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   9
  10 // Use, modification and distribution is subject to the Boost Software License,
  11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  12 // http://www.boost.org/LICENSE_1_0.txt)
  13
  14 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  15 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  16 // PROJ4 is maintained by Frank Warmerdam
  17 // PROJ4 is converted to Geometry Library by Barend Gehrels (Geodan, Amsterdam)
  18
  19 // Original copyright notice:
  20
  21 // Permission is hereby granted, free of charge, to any person obtaining a
  22 // copy of this software and associated documentation files (the "Software"),
  23 // to deal in the Software without restriction, including without limitation
  24 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  25 // and/or sell copies of the Software, and to permit persons to whom the
  26 // Software is furnished to do so, subject to the following conditions:
  27
  28 // The above copyright notice and this permission notice shall be included
  29 // in all copies or substantial portions of the Software.
  30
  31 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  32 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  33 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  34 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  35 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  36 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  37 // DEALINGS IN THE SOFTWARE.
  38
  39 #ifndef BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP
  40 #define BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP
  41
  42 #include <cstdlib>
  43 #include <string>
  44 #include <vector>
  45
  46 #include <boost/algorithm/string.hpp>
  47 #include <boost/lexical_cast.hpp>
  48 #include <boost/range.hpp>
  49 #include <boost/type_traits/is_same.hpp>
  50
  51 #include <boost/geometry/util/math.hpp>
  52 #include <boost/geometry/util/condition.hpp>
  53
  54 #include <boost/geometry/srs/projections/impl/dms_parser.hpp>
  55 #include <boost/geometry/srs/projections/impl/pj_datum_set.hpp>
  56 #include <boost/geometry/srs/projections/impl/pj_datums.hpp>
  57 #include <boost/geometry/srs/projections/impl/pj_ell_set.hpp>
  58 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  59 #include <boost/geometry/srs/projections/impl/pj_units.hpp>
  60 #include <boost/geometry/srs/projections/impl/projects.hpp>
  61 #include <boost/geometry/srs/projections/proj4.hpp>
  62
  63
  64 namespace boost { namespace geometry { namespace projections
  65 {
  66
  67
  68 namespace detail
  69 {
  70
  71 template <typename BGParams, typename T>
  72 inline void pj_push_defaults(BGParams const& /*bg_params*/, parameters<T>& pin)
  73 {
  74     pin.params.push_back(pj_mkparam<T>("ellps=WGS84"));
  75
  76     if (pin.name == "aea")
  77     {
  78         pin.params.push_back(pj_mkparam<T>("lat_1=29.5"));
  79         pin.params.push_back(pj_mkparam<T>("lat_2=45.5 "));
  80     }
  81     else if (pin.name == "lcc")
  82     {
  83         pin.params.push_back(pj_mkparam<T>("lat_1=33"));
  84         pin.params.push_back(pj_mkparam<T>("lat_2=45"));
  85     }
  86     else if (pin.name == "lagrng")
  87     {
  88         pin.params.push_back(pj_mkparam<T>("W=2"));
  89     }
  90 }
  91
  92 template <BOOST_GEOMETRY_PROJECTIONS_DETAIL_TYPENAME_PX, typename T>
  93 inline void pj_push_defaults(srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> const& /*bg_params*/,
  94                              parameters<T>& pin)
  95 {
  96     typedef srs::static_proj4<BOOST_GEOMETRY_PROJECTIONS_DETAIL_PX> static_parameters_type;
  97     typedef typename srs::par4::detail::pick_proj_tag
  98         <
  99             static_parameters_type
 100         >::type proj_tag;
 101
 102     // statically defaulting to WGS84
 103     //pin.params.push_back(pj_mkparam("ellps=WGS84"));
 104
 105     if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::aea>::value)))
 106     {
 107         pin.params.push_back(pj_mkparam<T>("lat_1=29.5"));
 108         pin.params.push_back(pj_mkparam<T>("lat_2=45.5 "));
 109     }
 110     else if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::lcc>::value)))
 111     {
 112         pin.params.push_back(pj_mkparam<T>("lat_1=33"));
 113         pin.params.push_back(pj_mkparam<T>("lat_2=45"));
 114     }
 115     else if (BOOST_GEOMETRY_CONDITION((boost::is_same<proj_tag, srs::par4::lagrng>::value)))
 116     {
 117         pin.params.push_back(pj_mkparam<T>("W=2"));
 118     }
 119 }
 120
 121 template <typename T>
 122 inline void pj_init_units(std::vector<pvalue<T> > const& params,
 123                           std::string const& sunits,
 124                           std::string const& sto_meter,
 125                           T & to_meter,
 126                           T & fr_meter,
 127                           T const& default_to_meter,
 128                           T const& default_fr_meter)
 129 {
 130     std::string s;
 131     std::string units = pj_param(params, sunits).s;
 132     if (! units.empty())
 133     {
 134         const int n = sizeof(pj_units) / sizeof(pj_units[0]);
 135         int index = -1;
 136         for (int i = 0; i < n && index == -1; i++)
 137         {
 138             if(pj_units[i].id == units)
 139             {
 140                 index = i;
 141             }
 142         }
 143
 144         if (index == -1) {
 145             BOOST_THROW_EXCEPTION( projection_exception(-7) );
 146         }
 147         s = pj_units[index].to_meter;
 148     }
 149
 150     if (s.empty())
 151     {
 152         s = pj_param(params, sto_meter).s;
 153     }
 154
 155     if (! s.empty())
 156     {
 157         std::size_t const pos = s.find('/');
 158         if (pos == std::string::npos)
 159         {
 160             to_meter = lexical_cast<T>(s);
 161         }
 162         else
 163         {
 164             T const numerator = lexical_cast<T>(s.substr(0, pos));
 165             T const denominator = lexical_cast<T>(s.substr(pos + 1));
 166             if (numerator == 0.0 || denominator == 0.0)
 167             {
 168                 BOOST_THROW_EXCEPTION( projection_exception(-99) );
 169             }
 170             to_meter = numerator / denominator;
 171         }
 172         if (to_meter == 0.0)
 173         {
 174             BOOST_THROW_EXCEPTION( projection_exception(-99) );
 175         }
 176         fr_meter = 1. / to_meter;
 177     }
 178     else
 179     {
 180         to_meter = default_to_meter;
 181         fr_meter = default_fr_meter;
 182     }
 183 }
 184
 185 /************************************************************************/
 186 /*                              pj_init()                               */
 187 /*                                                                      */
 188 /*      Main entry point for initialing a PJ projections                */
 189 /*      definition.  Note that the projection specific function is      */
 190 /*      called to do the initial allocation so it can be created        */
 191 /*      large enough to hold projection specific parameters.            */
 192 /************************************************************************/
 193 template <typename T, typename BGParams, typename R>
 194 inline parameters<T> pj_init(BGParams const& bg_params, R const& arguments, bool use_defaults = true)
 195 {
 196     parameters<T> pin;
 197     for (std::vector<std::string>::const_iterator it = boost::begin(arguments);
 198         it != boost::end(arguments); it++)
 199     {
 200         pin.params.push_back(pj_mkparam<T>(*it));
 201     }
 202
 203     /* check if +init present */
 204     if (pj_param(pin.params, "tinit").i)
 205     {
 206         // maybe TODO: handle "init" parameter
 207         //if (!(curr = get_init(&arguments, curr, pj_param(pin.params, "sinit").s)))
 208     }
 209
 210     // find projection -> implemented in projection factory
 211     pin.name = pj_param(pin.params, "sproj").s;
 212     // exception thrown in projection<>
 213     // TODO: consider throwing here both projection_unknown_id_exception and
 214     // projection_not_named_exception in order to throw before other exceptions
 215     //if (pin.name.empty())
 216     //{ BOOST_THROW_EXCEPTION( projection_not_named_exception() ); }
 217
 218     // set defaults, unless inhibited
 219     // GL-Addition, if use_defaults is false then defaults are ignored
 220     if (use_defaults && ! pj_param(pin.params, "bno_defs").i)
 221     {
 222         // proj4 gets defaults from "proj_def.dat", file of 94/02/23 with a few defaults.
 223         // Here manually
 224         pj_push_defaults(bg_params, pin);
 225         //curr = get_defaults(&arguments, curr, name);
 226     }
 227
 228     /* allocate projection structure */
 229     // done by BGParams constructor:
 230     // pin.is_latlong = 0;
 231     // pin.is_geocent = 0;
 232     // pin.long_wrap_center = 0.0;
 233     // pin.long_wrap_center = 0.0;
 234     pin.is_long_wrap_set = false;
 235
 236     /* set datum parameters */
 237     pj_datum_set(bg_params, pin.params, pin);
 238
 239     /* set ellipsoid/sphere parameters */
 240     pj_ell_set(bg_params, pin.params, pin.a, pin.es);
 241
 242     pin.a_orig = pin.a;
 243     pin.es_orig = pin.es;
 244
 245     pin.e = sqrt(pin.es);
 246     pin.ra = 1. / pin.a;
 247     pin.one_es = 1. - pin.es;
 248     if (pin.one_es == 0.) {
 249         BOOST_THROW_EXCEPTION( projection_exception(-6) );
 250     }
 251     pin.rone_es = 1./pin.one_es;
 252
 253     /* Now that we have ellipse information check for WGS84 datum */
 254     if( pin.datum_type == PJD_3PARAM
 255         && pin.datum_params[0] == 0.0
 256         && pin.datum_params[1] == 0.0
 257         && pin.datum_params[2] == 0.0
 258         && pin.a == 6378137.0
 259         && geometry::math::abs(pin.es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/
 260     {
 261         pin.datum_type = PJD_WGS84;
 262     }
 263
 264     /* set pin.geoc coordinate system */
 265     pin.geoc = (pin.es && pj_param(pin.params, "bgeoc").i);
 266
 267     /* over-ranging flag */
 268     pin.over = pj_param(pin.params, "bover").i;
 269
 270     /* longitude center for wrapping */
 271     pin.is_long_wrap_set = pj_param(pin.params, "tlon_wrap").i != 0;
 272     if (pin.is_long_wrap_set)
 273         pin.long_wrap_center = pj_param(pin.params, "rlon_wrap").f;
 274
 275     /* central meridian */
 276     pin.lam0 = pj_param(pin.params, "rlon_0").f;
 277
 278     /* central latitude */
 279     pin.phi0 = pj_param(pin.params, "rlat_0").f;
 280
 281     /* false easting and northing */
 282     pin.x0 = pj_param(pin.params, "dx_0").f;
 283     pin.y0 = pj_param(pin.params, "dy_0").f;
 284
 285     /* general scaling factor */
 286     if (pj_param(pin.params, "tk_0").i)
 287         pin.k0 = pj_param(pin.params, "dk_0").f;
 288     else if (pj_param(pin.params, "tk").i)
 289         pin.k0 = pj_param(pin.params, "dk").f;
 290     else
 291         pin.k0 = 1.;
 292     if (pin.k0 <= 0.) {
 293         BOOST_THROW_EXCEPTION( projection_exception(-31) );
 294     }
 295
 296     /* set units */
 297     pj_init_units(pin.params, "sunits", "sto_meter",
 298                   pin.to_meter, pin.fr_meter, 1., 1.);
 299     pj_init_units(pin.params, "svunits", "svto_meter",
 300                   pin.vto_meter, pin.vfr_meter, pin.to_meter, pin.fr_meter);
 301
 302     /* prime meridian */
 303     std::string pm = pj_param(pin.params, "spm").s;
 304     if (! pm.empty())
 305     {
 306         std::string value;
 307
 308         int n = sizeof(pj_prime_meridians) / sizeof(pj_prime_meridians[0]);
 309         for (int i = 0; i < n ; i++)
 310         {
 311             if(pj_prime_meridians[i].id == pm)
 312             {
 313                 value = pj_prime_meridians[i].defn;
 314                 break;
 315             }
 316         }
 317
 318         dms_parser<T, true> parser;
 319
 320         // TODO: Handle case when lexical_cast is not used consistently.
 321         //       This should probably be done in dms_parser.
 322         BOOST_TRY
 323         {
 324             if (value.empty()) {
 325                 pin.from_greenwich = parser.apply(pm).angle();
 326             } else {
 327                 pin.from_greenwich = parser.apply(value).angle();
 328             }
 329         }
 330         BOOST_CATCH(boost::bad_lexical_cast const&)
 331         {
 332             BOOST_THROW_EXCEPTION( projection_exception(-46) );
 333         }
 334         BOOST_CATCH_END
 335     }
 336     else
 337     {
 338         pin.from_greenwich = 0.0;
 339     }
 340
 341     return pin;
 342 }
 343
 344 /************************************************************************/
 345 /*                            pj_init_plus()                            */
 346 /*                                                                      */
 347 /*      Same as pj_init() except it takes one argument string with      */
 348 /*      individual arguments preceeded by '+', such as "+proj=utm       */
 349 /*      +zone=11 +ellps=WGS84".                                         */
 350 /************************************************************************/
 351 template <typename T, typename BGParams>
 352 inline parameters<T> pj_init_plus(BGParams const& bg_params, std::string const& definition, bool use_defaults = true)
 353 {
 354     const char* sep = " +";
 355
 356     /* split into arguments based on '+' and trim white space */
 357
 358     // boost::split splits on one character, here it should be on " +", so implementation below
 359     // todo: put in different routine or sort out
 360     std::vector<std::string> arguments;
 361     std::string def = boost::trim_copy(definition);
 362     boost::trim_left_if(def, boost::is_any_of(sep));
 363
 364     std::string::size_type loc = def.find(sep);
 365     while (loc != std::string::npos)
 366     {
 367         std::string par = def.substr(0, loc);
 368         boost::trim(par);
 369         if (! par.empty())
 370         {
 371             arguments.push_back(par);
 372         }
 373
 374         def.erase(0, loc);
 375         boost::trim_left_if(def, boost::is_any_of(sep));
 376         loc = def.find(sep);
 377     }
 378
 379     if (! def.empty())
 380     {
 381         arguments.push_back(def);
 382     }
 383
 384     /*boost::split(arguments, definition, boost::is_any_of("+"));
 385     for (std::vector<std::string>::iterator it = arguments.begin(); it != arguments.end(); it++)
 386     {
 387         boost::trim(*it);
 388     }*/
 389     return pj_init<T>(bg_params, arguments, use_defaults);
 390 }
 391
 392 } // namespace detail
 393 }}} // namespace boost::geometry::projections
 394
 395 #endif // BOOST_GEOMETRY_PROJECTIONS_IMPL_PJ_INIT_HPP