sconics.hpp 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568
  1. // Boost.Geometry - gis-projections (based on PROJ4)
  2. // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
  3. // This file was modified by Oracle on 2017, 2018.
  4. // Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
  5. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
  6. // Use, modification and distribution is subject to the Boost Software License,
  7. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8. // http://www.boost.org/LICENSE_1_0.txt)
  9. // This file is converted from PROJ4, http://trac.osgeo.org/proj
  10. // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  11. // PROJ4 is maintained by Frank Warmerdam
  12. // PROJ4 is converted to Boost.Geometry by Barend Gehrels
  13. // Last updated version of proj: 5.0.0
  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_SCONICS_HPP
  31. #define BOOST_GEOMETRY_PROJECTIONS_SCONICS_HPP
  32. #include <boost/geometry/util/math.hpp>
  33. #include <boost/math/special_functions/hypot.hpp>
  34. #include <boost/geometry/srs/projections/impl/base_static.hpp>
  35. #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  36. #include <boost/geometry/srs/projections/impl/projects.hpp>
  37. #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  38. namespace boost { namespace geometry
  39. {
  40. namespace srs { namespace par4
  41. {
  42. struct euler {}; // Euler
  43. struct murd1 {}; // Murdoch I
  44. struct murd2 {}; // Murdoch II
  45. struct murd3 {}; // Murdoch III
  46. struct pconic {}; // Perspective Conic
  47. struct tissot {}; // Tissot
  48. struct vitk1 {}; // Vitkovsky I
  49. }} //namespace srs::par4
  50. namespace projections
  51. {
  52. #ifndef DOXYGEN_NO_DETAIL
  53. namespace detail { namespace sconics
  54. {
  55. enum proj_type {
  56. proj_euler = 0,
  57. proj_murd1 = 1,
  58. proj_murd2 = 2,
  59. proj_murd3 = 3,
  60. proj_pconic = 4,
  61. proj_tissot = 5,
  62. proj_vitk1 = 6
  63. };
  64. static const double epsilon10 = 1.e-10;
  65. static const double epsilon = 1e-10;
  66. template <typename T>
  67. struct par_sconics
  68. {
  69. T n;
  70. T rho_c;
  71. T rho_0;
  72. T sig;
  73. T c1, c2;
  74. proj_type type;
  75. };
  76. /* get common factors for simple conics */
  77. template <typename Parameters, typename T>
  78. inline int phi12(Parameters& par, par_sconics<T>& proj_parm, T *del)
  79. {
  80. T p1, p2;
  81. int err = 0;
  82. if (!pj_param_r(par.params, "lat_1", p1) ||
  83. !pj_param_r(par.params, "lat_2", p2)) {
  84. err = -41;
  85. } else {
  86. //p1 = pj_get_param_r(par.params, "lat_1"); // set above
  87. //p2 = pj_get_param_r(par.params, "lat_2"); // set above
  88. *del = 0.5 * (p2 - p1);
  89. proj_parm.sig = 0.5 * (p2 + p1);
  90. err = (fabs(*del) < epsilon || fabs(proj_parm.sig) < epsilon) ? -42 : 0;
  91. }
  92. return err;
  93. }
  94. // template class, using CRTP to implement forward/inverse
  95. template <typename T, typename Parameters>
  96. struct base_sconics_spheroid
  97. : public base_t_fi<base_sconics_spheroid<T, Parameters>, T, Parameters>
  98. {
  99. par_sconics<T> m_proj_parm;
  100. inline base_sconics_spheroid(const Parameters& par)
  101. : base_t_fi<base_sconics_spheroid<T, Parameters>, T, Parameters>(*this, par)
  102. {}
  103. // FORWARD(s_forward) spheroid
  104. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  105. inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const
  106. {
  107. T rho;
  108. switch (this->m_proj_parm.type) {
  109. case proj_murd2:
  110. rho = this->m_proj_parm.rho_c + tan(this->m_proj_parm.sig - lp_lat);
  111. break;
  112. case proj_pconic:
  113. rho = this->m_proj_parm.c2 * (this->m_proj_parm.c1 - tan(lp_lat - this->m_proj_parm.sig));
  114. break;
  115. default:
  116. rho = this->m_proj_parm.rho_c - lp_lat;
  117. break;
  118. }
  119. xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
  120. xy_y = this->m_proj_parm.rho_0 - rho * cos(lp_lon);
  121. }
  122. // INVERSE(s_inverse) ellipsoid & spheroid
  123. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  124. inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const
  125. {
  126. T rho;
  127. rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho_0 - xy_y);
  128. if (this->m_proj_parm.n < 0.) {
  129. rho = - rho;
  130. xy_x = - xy_x;
  131. xy_y = - xy_y;
  132. }
  133. lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
  134. switch (this->m_proj_parm.type) {
  135. case proj_pconic:
  136. lp_lat = atan(this->m_proj_parm.c1 - rho / this->m_proj_parm.c2) + this->m_proj_parm.sig;
  137. break;
  138. case proj_murd2:
  139. lp_lat = this->m_proj_parm.sig - atan(rho - this->m_proj_parm.rho_c);
  140. break;
  141. default:
  142. lp_lat = this->m_proj_parm.rho_c - rho;
  143. }
  144. }
  145. static inline std::string get_name()
  146. {
  147. return "sconics_spheroid";
  148. }
  149. };
  150. template <typename Parameters, typename T>
  151. inline void setup(Parameters& par, par_sconics<T>& proj_parm, proj_type type)
  152. {
  153. static const T half_pi = detail::half_pi<T>();
  154. T del, cs;
  155. int err;
  156. proj_parm.type = type;
  157. err = phi12(par, proj_parm, &del);
  158. if(err)
  159. BOOST_THROW_EXCEPTION( projection_exception(err) );
  160. switch (proj_parm.type) {
  161. case proj_tissot:
  162. proj_parm.n = sin(proj_parm.sig);
  163. cs = cos(del);
  164. proj_parm.rho_c = proj_parm.n / cs + cs / proj_parm.n;
  165. proj_parm.rho_0 = sqrt((proj_parm.rho_c - 2 * sin(par.phi0))/proj_parm.n);
  166. break;
  167. case proj_murd1:
  168. proj_parm.rho_c = sin(del)/(del * tan(proj_parm.sig)) + proj_parm.sig;
  169. proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
  170. proj_parm.n = sin(proj_parm.sig);
  171. break;
  172. case proj_murd2:
  173. proj_parm.rho_c = (cs = sqrt(cos(del))) / tan(proj_parm.sig);
  174. proj_parm.rho_0 = proj_parm.rho_c + tan(proj_parm.sig - par.phi0);
  175. proj_parm.n = sin(proj_parm.sig) * cs;
  176. break;
  177. case proj_murd3:
  178. proj_parm.rho_c = del / (tan(proj_parm.sig) * tan(del)) + proj_parm.sig;
  179. proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
  180. proj_parm.n = sin(proj_parm.sig) * sin(del) * tan(del) / (del * del);
  181. break;
  182. case proj_euler:
  183. proj_parm.n = sin(proj_parm.sig) * sin(del) / del;
  184. del *= 0.5;
  185. proj_parm.rho_c = del / (tan(del) * tan(proj_parm.sig)) + proj_parm.sig;
  186. proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
  187. break;
  188. case proj_pconic:
  189. proj_parm.n = sin(proj_parm.sig);
  190. proj_parm.c2 = cos(del);
  191. proj_parm.c1 = 1./tan(proj_parm.sig);
  192. if (fabs(del = par.phi0 - proj_parm.sig) - epsilon10 >= half_pi)
  193. BOOST_THROW_EXCEPTION( projection_exception(error_lat_0_half_pi_from_mean) );
  194. proj_parm.rho_0 = proj_parm.c2 * (proj_parm.c1 - tan(del));
  195. break;
  196. case proj_vitk1:
  197. proj_parm.n = (cs = tan(del)) * sin(proj_parm.sig) / del;
  198. proj_parm.rho_c = del / (cs * tan(proj_parm.sig)) + proj_parm.sig;
  199. proj_parm.rho_0 = proj_parm.rho_c - par.phi0;
  200. break;
  201. }
  202. par.es = 0;
  203. }
  204. // Euler
  205. template <typename Parameters, typename T>
  206. inline void setup_euler(Parameters& par, par_sconics<T>& proj_parm)
  207. {
  208. setup(par, proj_parm, proj_euler);
  209. }
  210. // Tissot
  211. template <typename Parameters, typename T>
  212. inline void setup_tissot(Parameters& par, par_sconics<T>& proj_parm)
  213. {
  214. setup(par, proj_parm, proj_tissot);
  215. }
  216. // Murdoch I
  217. template <typename Parameters, typename T>
  218. inline void setup_murd1(Parameters& par, par_sconics<T>& proj_parm)
  219. {
  220. setup(par, proj_parm, proj_murd1);
  221. }
  222. // Murdoch II
  223. template <typename Parameters, typename T>
  224. inline void setup_murd2(Parameters& par, par_sconics<T>& proj_parm)
  225. {
  226. setup(par, proj_parm, proj_murd2);
  227. }
  228. // Murdoch III
  229. template <typename Parameters, typename T>
  230. inline void setup_murd3(Parameters& par, par_sconics<T>& proj_parm)
  231. {
  232. setup(par, proj_parm, proj_murd3);
  233. }
  234. // Perspective Conic
  235. template <typename Parameters, typename T>
  236. inline void setup_pconic(Parameters& par, par_sconics<T>& proj_parm)
  237. {
  238. setup(par, proj_parm, proj_pconic);
  239. }
  240. // Vitkovsky I
  241. template <typename Parameters, typename T>
  242. inline void setup_vitk1(Parameters& par, par_sconics<T>& proj_parm)
  243. {
  244. setup(par, proj_parm, proj_vitk1);
  245. }
  246. }} // namespace detail::sconics
  247. #endif // doxygen
  248. /*!
  249. \brief Tissot projection
  250. \ingroup projections
  251. \tparam Geographic latlong point type
  252. \tparam Cartesian xy point type
  253. \tparam Parameters parameter type
  254. \par Projection characteristics
  255. - Conic
  256. - Spheroid
  257. \par Projection parameters
  258. - lat_1: Latitude of first standard parallel
  259. - lat_2: Latitude of second standard parallel
  260. \par Example
  261. \image html ex_tissot.gif
  262. */
  263. template <typename T, typename Parameters>
  264. struct tissot_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
  265. {
  266. inline tissot_spheroid(const Parameters& par) : detail::sconics::base_sconics_spheroid<T, Parameters>(par)
  267. {
  268. detail::sconics::setup_tissot(this->m_par, this->m_proj_parm);
  269. }
  270. };
  271. /*!
  272. \brief Murdoch I projection
  273. \ingroup projections
  274. \tparam Geographic latlong point type
  275. \tparam Cartesian xy point type
  276. \tparam Parameters parameter type
  277. \par Projection characteristics
  278. - Conic
  279. - Spheroid
  280. \par Projection parameters
  281. - lat_1: Latitude of first standard parallel
  282. - lat_2: Latitude of second standard parallel
  283. \par Example
  284. \image html ex_murd1.gif
  285. */
  286. template <typename T, typename Parameters>
  287. struct murd1_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
  288. {
  289. inline murd1_spheroid(const Parameters& par) : detail::sconics::base_sconics_spheroid<T, Parameters>(par)
  290. {
  291. detail::sconics::setup_murd1(this->m_par, this->m_proj_parm);
  292. }
  293. };
  294. /*!
  295. \brief Murdoch II projection
  296. \ingroup projections
  297. \tparam Geographic latlong point type
  298. \tparam Cartesian xy point type
  299. \tparam Parameters parameter type
  300. \par Projection characteristics
  301. - Conic
  302. - Spheroid
  303. \par Projection parameters
  304. - lat_1: Latitude of first standard parallel
  305. - lat_2: Latitude of second standard parallel
  306. \par Example
  307. \image html ex_murd2.gif
  308. */
  309. template <typename T, typename Parameters>
  310. struct murd2_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
  311. {
  312. inline murd2_spheroid(const Parameters& par) : detail::sconics::base_sconics_spheroid<T, Parameters>(par)
  313. {
  314. detail::sconics::setup_murd2(this->m_par, this->m_proj_parm);
  315. }
  316. };
  317. /*!
  318. \brief Murdoch III projection
  319. \ingroup projections
  320. \tparam Geographic latlong point type
  321. \tparam Cartesian xy point type
  322. \tparam Parameters parameter type
  323. \par Projection characteristics
  324. - Conic
  325. - Spheroid
  326. \par Projection parameters
  327. - lat_1: Latitude of first standard parallel
  328. - lat_2: Latitude of second standard parallel
  329. \par Example
  330. \image html ex_murd3.gif
  331. */
  332. template <typename T, typename Parameters>
  333. struct murd3_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
  334. {
  335. inline murd3_spheroid(const Parameters& par) : detail::sconics::base_sconics_spheroid<T, Parameters>(par)
  336. {
  337. detail::sconics::setup_murd3(this->m_par, this->m_proj_parm);
  338. }
  339. };
  340. /*!
  341. \brief Euler projection
  342. \ingroup projections
  343. \tparam Geographic latlong point type
  344. \tparam Cartesian xy point type
  345. \tparam Parameters parameter type
  346. \par Projection characteristics
  347. - Conic
  348. - Spheroid
  349. \par Projection parameters
  350. - lat_1: Latitude of first standard parallel
  351. - lat_2: Latitude of second standard parallel
  352. \par Example
  353. \image html ex_euler.gif
  354. */
  355. template <typename T, typename Parameters>
  356. struct euler_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
  357. {
  358. inline euler_spheroid(const Parameters& par) : detail::sconics::base_sconics_spheroid<T, Parameters>(par)
  359. {
  360. detail::sconics::setup_euler(this->m_par, this->m_proj_parm);
  361. }
  362. };
  363. /*!
  364. \brief Perspective Conic projection
  365. \ingroup projections
  366. \tparam Geographic latlong point type
  367. \tparam Cartesian xy point type
  368. \tparam Parameters parameter type
  369. \par Projection characteristics
  370. - Conic
  371. - Spheroid
  372. \par Projection parameters
  373. - lat_1: Latitude of first standard parallel
  374. - lat_2: Latitude of second standard parallel
  375. \par Example
  376. \image html ex_pconic.gif
  377. */
  378. template <typename T, typename Parameters>
  379. struct pconic_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
  380. {
  381. inline pconic_spheroid(const Parameters& par) : detail::sconics::base_sconics_spheroid<T, Parameters>(par)
  382. {
  383. detail::sconics::setup_pconic(this->m_par, this->m_proj_parm);
  384. }
  385. };
  386. /*!
  387. \brief Vitkovsky I projection
  388. \ingroup projections
  389. \tparam Geographic latlong point type
  390. \tparam Cartesian xy point type
  391. \tparam Parameters parameter type
  392. \par Projection characteristics
  393. - Conic
  394. - Spheroid
  395. \par Projection parameters
  396. - lat_1: Latitude of first standard parallel
  397. - lat_2: Latitude of second standard parallel
  398. \par Example
  399. \image html ex_vitk1.gif
  400. */
  401. template <typename T, typename Parameters>
  402. struct vitk1_spheroid : public detail::sconics::base_sconics_spheroid<T, Parameters>
  403. {
  404. inline vitk1_spheroid(const Parameters& par) : detail::sconics::base_sconics_spheroid<T, Parameters>(par)
  405. {
  406. detail::sconics::setup_vitk1(this->m_par, this->m_proj_parm);
  407. }
  408. };
  409. #ifndef DOXYGEN_NO_DETAIL
  410. namespace detail
  411. {
  412. // Static projection
  413. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::euler, euler_spheroid, euler_spheroid)
  414. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::murd1, murd1_spheroid, murd1_spheroid)
  415. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::murd2, murd2_spheroid, murd2_spheroid)
  416. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::murd3, murd3_spheroid, murd3_spheroid)
  417. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::pconic, pconic_spheroid, pconic_spheroid)
  418. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::tissot, tissot_spheroid, tissot_spheroid)
  419. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::vitk1, vitk1_spheroid, vitk1_spheroid)
  420. // Factory entry(s)
  421. template <typename T, typename Parameters>
  422. class tissot_entry : public detail::factory_entry<T, Parameters>
  423. {
  424. public :
  425. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  426. {
  427. return new base_v_fi<tissot_spheroid<T, Parameters>, T, Parameters>(par);
  428. }
  429. };
  430. template <typename T, typename Parameters>
  431. class murd1_entry : public detail::factory_entry<T, Parameters>
  432. {
  433. public :
  434. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  435. {
  436. return new base_v_fi<murd1_spheroid<T, Parameters>, T, Parameters>(par);
  437. }
  438. };
  439. template <typename T, typename Parameters>
  440. class murd2_entry : public detail::factory_entry<T, Parameters>
  441. {
  442. public :
  443. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  444. {
  445. return new base_v_fi<murd2_spheroid<T, Parameters>, T, Parameters>(par);
  446. }
  447. };
  448. template <typename T, typename Parameters>
  449. class murd3_entry : public detail::factory_entry<T, Parameters>
  450. {
  451. public :
  452. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  453. {
  454. return new base_v_fi<murd3_spheroid<T, Parameters>, T, Parameters>(par);
  455. }
  456. };
  457. template <typename T, typename Parameters>
  458. class euler_entry : public detail::factory_entry<T, Parameters>
  459. {
  460. public :
  461. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  462. {
  463. return new base_v_fi<euler_spheroid<T, Parameters>, T, Parameters>(par);
  464. }
  465. };
  466. template <typename T, typename Parameters>
  467. class pconic_entry : public detail::factory_entry<T, Parameters>
  468. {
  469. public :
  470. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  471. {
  472. return new base_v_fi<pconic_spheroid<T, Parameters>, T, Parameters>(par);
  473. }
  474. };
  475. template <typename T, typename Parameters>
  476. class vitk1_entry : public detail::factory_entry<T, Parameters>
  477. {
  478. public :
  479. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  480. {
  481. return new base_v_fi<vitk1_spheroid<T, Parameters>, T, Parameters>(par);
  482. }
  483. };
  484. template <typename T, typename Parameters>
  485. inline void sconics_init(detail::base_factory<T, Parameters>& factory)
  486. {
  487. factory.add_to_factory("tissot", new tissot_entry<T, Parameters>);
  488. factory.add_to_factory("murd1", new murd1_entry<T, Parameters>);
  489. factory.add_to_factory("murd2", new murd2_entry<T, Parameters>);
  490. factory.add_to_factory("murd3", new murd3_entry<T, Parameters>);
  491. factory.add_to_factory("euler", new euler_entry<T, Parameters>);
  492. factory.add_to_factory("pconic", new pconic_entry<T, Parameters>);
  493. factory.add_to_factory("vitk1", new vitk1_entry<T, Parameters>);
  494. }
  495. } // namespace detail
  496. #endif // doxygen
  497. } // namespace projections
  498. }} // namespace boost::geometry
  499. #endif // BOOST_GEOMETRY_PROJECTIONS_SCONICS_HPP