pj_init.hpp 13 KB

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