sconics.hpp 19 KB

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