gn_sinu.hpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423
  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_GN_SINU_HPP
  31. #define BOOST_GEOMETRY_PROJECTIONS_GN_SINU_HPP
  32. #include <boost/geometry/util/math.hpp>
  33. #include <boost/geometry/srs/projections/impl/base_static.hpp>
  34. #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  35. #include <boost/geometry/srs/projections/impl/projects.hpp>
  36. #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  37. #include <boost/geometry/srs/projections/impl/aasincos.hpp>
  38. #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
  39. namespace boost { namespace geometry
  40. {
  41. namespace srs { namespace par4
  42. {
  43. struct gn_sinu {}; // General Sinusoidal Series
  44. struct sinu {}; // Sinusoidal (Sanson-Flamsteed)
  45. struct eck6 {}; // Eckert VI
  46. struct mbtfps {}; // McBryde-Thomas Flat-Polar Sinusoidal
  47. }} //namespace srs::par4
  48. namespace projections
  49. {
  50. #ifndef DOXYGEN_NO_DETAIL
  51. namespace detail { namespace gn_sinu
  52. {
  53. static const double epsilon10 = 1e-10;
  54. static const int max_iter = 8;
  55. static const double loop_tol = 1e-7;
  56. template <typename T>
  57. struct par_gn_sinu
  58. {
  59. detail::en<T> en;
  60. T m, n, C_x, C_y;
  61. };
  62. /* Ellipsoidal Sinusoidal only */
  63. // template class, using CRTP to implement forward/inverse
  64. template <typename T, typename Parameters>
  65. struct base_gn_sinu_ellipsoid
  66. : public base_t_fi<base_gn_sinu_ellipsoid<T, Parameters>, T, Parameters>
  67. {
  68. par_gn_sinu<T> m_proj_parm;
  69. inline base_gn_sinu_ellipsoid(const Parameters& par)
  70. : base_t_fi<base_gn_sinu_ellipsoid<T, Parameters>, T, Parameters>(*this, par)
  71. {}
  72. // FORWARD(e_forward) ellipsoid
  73. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  74. inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const
  75. {
  76. T s, c;
  77. xy_y = pj_mlfn(lp_lat, s = sin(lp_lat), c = cos(lp_lat), this->m_proj_parm.en);
  78. xy_x = lp_lon * c / sqrt(1. - this->m_par.es * s * s);
  79. }
  80. // INVERSE(e_inverse) ellipsoid
  81. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  82. inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const
  83. {
  84. static const T half_pi = detail::half_pi<T>();
  85. T s;
  86. if ((s = fabs(lp_lat = pj_inv_mlfn(xy_y, this->m_par.es, this->m_proj_parm.en))) < half_pi) {
  87. s = sin(lp_lat);
  88. lp_lon = xy_x * sqrt(1. - this->m_par.es * s * s) / cos(lp_lat);
  89. } else if ((s - epsilon10) < half_pi)
  90. lp_lon = 0.;
  91. else
  92. BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
  93. }
  94. /* General spherical sinusoidals */
  95. static inline std::string get_name()
  96. {
  97. return "gn_sinu_ellipsoid";
  98. }
  99. };
  100. // template class, using CRTP to implement forward/inverse
  101. template <typename T, typename Parameters>
  102. struct base_gn_sinu_spheroid
  103. : public base_t_fi<base_gn_sinu_spheroid<T, Parameters>, T, Parameters>
  104. {
  105. par_gn_sinu<T> m_proj_parm;
  106. inline base_gn_sinu_spheroid(const Parameters& par)
  107. : base_t_fi<base_gn_sinu_spheroid<T, Parameters>, T, Parameters>(*this, par)
  108. {}
  109. // FORWARD(s_forward) sphere
  110. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  111. inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const
  112. {
  113. if (this->m_proj_parm.m == 0.0)
  114. lp_lat = this->m_proj_parm.n != 1. ? aasin(this->m_proj_parm.n * sin(lp_lat)): lp_lat;
  115. else {
  116. T k, V;
  117. int i;
  118. k = this->m_proj_parm.n * sin(lp_lat);
  119. for (i = max_iter; i ; --i) {
  120. lp_lat -= V = (this->m_proj_parm.m * lp_lat + sin(lp_lat) - k) /
  121. (this->m_proj_parm.m + cos(lp_lat));
  122. if (fabs(V) < loop_tol)
  123. break;
  124. }
  125. if (!i) {
  126. BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
  127. }
  128. }
  129. xy_x = this->m_proj_parm.C_x * lp_lon * (this->m_proj_parm.m + cos(lp_lat));
  130. xy_y = this->m_proj_parm.C_y * lp_lat;
  131. }
  132. // INVERSE(s_inverse) sphere
  133. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  134. inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const
  135. {
  136. xy_y /= this->m_proj_parm.C_y;
  137. lp_lat = (this->m_proj_parm.m != 0.0) ? aasin((this->m_proj_parm.m * xy_y + sin(xy_y)) / this->m_proj_parm.n) :
  138. ( this->m_proj_parm.n != 1. ? aasin(sin(xy_y) / this->m_proj_parm.n) : xy_y );
  139. lp_lon = xy_x / (this->m_proj_parm.C_x * (this->m_proj_parm.m + cos(xy_y)));
  140. }
  141. static inline std::string get_name()
  142. {
  143. return "gn_sinu_spheroid";
  144. }
  145. };
  146. template <typename Parameters, typename T>
  147. inline void setup(Parameters& par, par_gn_sinu<T>& proj_parm)
  148. {
  149. par.es = 0;
  150. proj_parm.C_x = (proj_parm.C_y = sqrt((proj_parm.m + 1.) / proj_parm.n))/(proj_parm.m + 1.);
  151. }
  152. // General Sinusoidal Series
  153. template <typename Parameters, typename T>
  154. inline void setup_gn_sinu(Parameters& par, par_gn_sinu<T>& proj_parm)
  155. {
  156. if (pj_param_f(par.params, "n", proj_parm.n)
  157. && pj_param_f(par.params, "m", proj_parm.m)) {
  158. if (proj_parm.n <= 0 || proj_parm.m < 0)
  159. BOOST_THROW_EXCEPTION( projection_exception(error_invalid_m_or_n) );
  160. } else
  161. BOOST_THROW_EXCEPTION( projection_exception(error_invalid_m_or_n) );
  162. setup(par, proj_parm);
  163. }
  164. // Sinusoidal (Sanson-Flamsteed)
  165. template <typename Parameters, typename T>
  166. inline void setup_sinu(Parameters& par, par_gn_sinu<T>& proj_parm)
  167. {
  168. proj_parm.en = pj_enfn<T>(par.es);
  169. if (par.es != 0.0) {
  170. /* empty */
  171. } else {
  172. proj_parm.n = 1.;
  173. proj_parm.m = 0.;
  174. setup(par, proj_parm);
  175. }
  176. }
  177. // Eckert VI
  178. template <typename Parameters, typename T>
  179. inline void setup_eck6(Parameters& par, par_gn_sinu<T>& proj_parm)
  180. {
  181. proj_parm.m = 1.;
  182. proj_parm.n = 2.570796326794896619231321691;
  183. setup(par, proj_parm);
  184. }
  185. // McBryde-Thomas Flat-Polar Sinusoidal
  186. template <typename Parameters, typename T>
  187. inline void setup_mbtfps(Parameters& par, par_gn_sinu<T>& proj_parm)
  188. {
  189. proj_parm.m = 0.5;
  190. proj_parm.n = 1.785398163397448309615660845;
  191. setup(par, proj_parm);
  192. }
  193. }} // namespace detail::gn_sinu
  194. #endif // doxygen
  195. /*!
  196. \brief General Sinusoidal Series projection
  197. \ingroup projections
  198. \tparam Geographic latlong point type
  199. \tparam Cartesian xy point type
  200. \tparam Parameters parameter type
  201. \par Projection characteristics
  202. - Pseudocylindrical
  203. - Spheroid
  204. \par Projection parameters
  205. - m (real)
  206. - n (real)
  207. \par Example
  208. \image html ex_gn_sinu.gif
  209. */
  210. template <typename T, typename Parameters>
  211. struct gn_sinu_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
  212. {
  213. inline gn_sinu_spheroid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>(par)
  214. {
  215. detail::gn_sinu::setup_gn_sinu(this->m_par, this->m_proj_parm);
  216. }
  217. };
  218. /*!
  219. \brief Sinusoidal (Sanson-Flamsteed) projection
  220. \ingroup projections
  221. \tparam Geographic latlong point type
  222. \tparam Cartesian xy point type
  223. \tparam Parameters parameter type
  224. \par Projection characteristics
  225. - Pseudocylindrical
  226. - Spheroid
  227. - Ellipsoid
  228. \par Example
  229. \image html ex_sinu.gif
  230. */
  231. template <typename T, typename Parameters>
  232. struct sinu_ellipsoid : public detail::gn_sinu::base_gn_sinu_ellipsoid<T, Parameters>
  233. {
  234. inline sinu_ellipsoid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_ellipsoid<T, Parameters>(par)
  235. {
  236. detail::gn_sinu::setup_sinu(this->m_par, this->m_proj_parm);
  237. }
  238. };
  239. /*!
  240. \brief Sinusoidal (Sanson-Flamsteed) 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. - Pseudocylindrical
  247. - Spheroid
  248. - Ellipsoid
  249. \par Example
  250. \image html ex_sinu.gif
  251. */
  252. template <typename T, typename Parameters>
  253. struct sinu_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
  254. {
  255. inline sinu_spheroid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>(par)
  256. {
  257. detail::gn_sinu::setup_sinu(this->m_par, this->m_proj_parm);
  258. }
  259. };
  260. /*!
  261. \brief Eckert VI projection
  262. \ingroup projections
  263. \tparam Geographic latlong point type
  264. \tparam Cartesian xy point type
  265. \tparam Parameters parameter type
  266. \par Projection characteristics
  267. - Pseudocylindrical
  268. - Spheroid
  269. \par Example
  270. \image html ex_eck6.gif
  271. */
  272. template <typename T, typename Parameters>
  273. struct eck6_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
  274. {
  275. inline eck6_spheroid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>(par)
  276. {
  277. detail::gn_sinu::setup_eck6(this->m_par, this->m_proj_parm);
  278. }
  279. };
  280. /*!
  281. \brief McBryde-Thomas Flat-Polar Sinusoidal projection
  282. \ingroup projections
  283. \tparam Geographic latlong point type
  284. \tparam Cartesian xy point type
  285. \tparam Parameters parameter type
  286. \par Projection characteristics
  287. - Pseudocylindrical
  288. - Spheroid
  289. \par Example
  290. \image html ex_mbtfps.gif
  291. */
  292. template <typename T, typename Parameters>
  293. struct mbtfps_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>
  294. {
  295. inline mbtfps_spheroid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_spheroid<T, Parameters>(par)
  296. {
  297. detail::gn_sinu::setup_mbtfps(this->m_par, this->m_proj_parm);
  298. }
  299. };
  300. #ifndef DOXYGEN_NO_DETAIL
  301. namespace detail
  302. {
  303. // Static projection
  304. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::gn_sinu, gn_sinu_spheroid, gn_sinu_spheroid)
  305. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::sinu, sinu_spheroid, sinu_ellipsoid)
  306. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::eck6, eck6_spheroid, eck6_spheroid)
  307. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::mbtfps, mbtfps_spheroid, mbtfps_spheroid)
  308. // Factory entry(s)
  309. template <typename T, typename Parameters>
  310. class gn_sinu_entry : public detail::factory_entry<T, Parameters>
  311. {
  312. public :
  313. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  314. {
  315. return new base_v_fi<gn_sinu_spheroid<T, Parameters>, T, Parameters>(par);
  316. }
  317. };
  318. template <typename T, typename Parameters>
  319. class sinu_entry : public detail::factory_entry<T, Parameters>
  320. {
  321. public :
  322. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  323. {
  324. if (par.es)
  325. return new base_v_fi<sinu_ellipsoid<T, Parameters>, T, Parameters>(par);
  326. else
  327. return new base_v_fi<sinu_spheroid<T, Parameters>, T, Parameters>(par);
  328. }
  329. };
  330. template <typename T, typename Parameters>
  331. class eck6_entry : public detail::factory_entry<T, Parameters>
  332. {
  333. public :
  334. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  335. {
  336. return new base_v_fi<eck6_spheroid<T, Parameters>, T, Parameters>(par);
  337. }
  338. };
  339. template <typename T, typename Parameters>
  340. class mbtfps_entry : public detail::factory_entry<T, Parameters>
  341. {
  342. public :
  343. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  344. {
  345. return new base_v_fi<mbtfps_spheroid<T, Parameters>, T, Parameters>(par);
  346. }
  347. };
  348. template <typename T, typename Parameters>
  349. inline void gn_sinu_init(detail::base_factory<T, Parameters>& factory)
  350. {
  351. factory.add_to_factory("gn_sinu", new gn_sinu_entry<T, Parameters>);
  352. factory.add_to_factory("sinu", new sinu_entry<T, Parameters>);
  353. factory.add_to_factory("eck6", new eck6_entry<T, Parameters>);
  354. factory.add_to_factory("mbtfps", new mbtfps_entry<T, Parameters>);
  355. }
  356. } // namespace detail
  357. #endif // doxygen
  358. } // namespace projections
  359. }} // namespace boost::geometry
  360. #endif // BOOST_GEOMETRY_PROJECTIONS_GN_SINU_HPP