gn_sinu.hpp 15 KB

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