sconics.hpp 20 KB

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