cass.hpp 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284
  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_CASS_HPP
  31. #define BOOST_GEOMETRY_PROJECTIONS_CASS_HPP
  32. #include <boost/geometry/srs/projections/impl/base_static.hpp>
  33. #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  34. #include <boost/geometry/srs/projections/impl/projects.hpp>
  35. #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  36. #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
  37. namespace boost { namespace geometry
  38. {
  39. namespace srs { namespace par4
  40. {
  41. struct cass {};
  42. }} //namespace srs::par4
  43. namespace projections
  44. {
  45. #ifndef DOXYGEN_NO_DETAIL
  46. namespace detail { namespace cass
  47. {
  48. //static const double epsilon10 = 1e-10;
  49. //static const double C1 = .16666666666666666666;
  50. //static const double C2 = .00833333333333333333;
  51. //static const double C3 = .04166666666666666666;
  52. //static const double C4 = .33333333333333333333;
  53. //static const double C5 = .06666666666666666666;
  54. template <typename T>
  55. inline T C1() { return .16666666666666666666666666666666666666; }
  56. template <typename T>
  57. inline T C2() { return .00833333333333333333333333333333333333; }
  58. template <typename T>
  59. inline T C3() { return .04166666666666666666666666666666666666; }
  60. template <typename T>
  61. inline T C4() { return .33333333333333333333333333333333333333; }
  62. template <typename T>
  63. inline T C5() { return .06666666666666666666666666666666666666; }
  64. template <typename T>
  65. struct par_cass
  66. {
  67. T m0;
  68. detail::en<T> en;
  69. };
  70. // template class, using CRTP to implement forward/inverse
  71. template <typename T, typename Parameters>
  72. struct base_cass_ellipsoid
  73. : public base_t_fi<base_cass_ellipsoid<T, Parameters>, T, Parameters>
  74. {
  75. par_cass<T> m_proj_parm;
  76. inline base_cass_ellipsoid(const Parameters& par)
  77. : base_t_fi<base_cass_ellipsoid<T, Parameters>, T, Parameters>(*this, par)
  78. {}
  79. // FORWARD(e_forward) ellipsoid
  80. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  81. inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const
  82. {
  83. static const T C1 = cass::C1<T>();
  84. static const T C2 = cass::C2<T>();
  85. static const T C3 = cass::C3<T>();
  86. T n = sin(lp_lat);
  87. T c = cos(lp_lat);
  88. xy_y = pj_mlfn(lp_lat, n, c, this->m_proj_parm.en);
  89. n = 1./sqrt(1. - this->m_par.es * n * n);
  90. T tn = tan(lp_lat);
  91. T t = tn * tn;
  92. T a1 = lp_lon * c;
  93. c *= this->m_par.es * c / (1 - this->m_par.es);
  94. T a2 = a1 * a1;
  95. xy_x = n * a1 * (1. - a2 * t *
  96. (C1 - (8. - t + 8. * c) * a2 * C2));
  97. xy_y -= this->m_proj_parm.m0 - n * tn * a2 *
  98. (.5 + (5. - t + 6. * c) * a2 * C3);
  99. }
  100. // INVERSE(e_inverse) ellipsoid
  101. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  102. inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const
  103. {
  104. static const T C3 = cass::C3<T>();
  105. static const T C4 = cass::C4<T>();
  106. static const T C5 = cass::C5<T>();
  107. T ph1;
  108. ph1 = pj_inv_mlfn(this->m_proj_parm.m0 + xy_y, this->m_par.es, this->m_proj_parm.en);
  109. T tn = tan(ph1); T t = tn * tn;
  110. T n = sin(ph1);
  111. T r = 1. / (1. - this->m_par.es * n * n);
  112. n = sqrt(r);
  113. r *= (1. - this->m_par.es) * n;
  114. T dd = xy_x / n;
  115. T d2 = dd * dd;
  116. lp_lat = ph1 - (n * tn / r) * d2 *
  117. (.5 - (1. + 3. * t) * d2 * C3);
  118. lp_lon = dd * (1. + t * d2 *
  119. (-C4 + (1. + 3. * t) * d2 * C5)) / cos(ph1);
  120. }
  121. static inline std::string get_name()
  122. {
  123. return "cass_ellipsoid";
  124. }
  125. };
  126. // template class, using CRTP to implement forward/inverse
  127. template <typename T, typename Parameters>
  128. struct base_cass_spheroid
  129. : public base_t_fi<base_cass_spheroid<T, Parameters>, T, Parameters>
  130. {
  131. par_cass<T> m_proj_parm;
  132. inline base_cass_spheroid(const Parameters& par)
  133. : base_t_fi<base_cass_spheroid<T, Parameters>, T, Parameters>(*this, par)
  134. {}
  135. // FORWARD(s_forward) spheroid
  136. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  137. inline void fwd(T& lp_lon, T& lp_lat, T& xy_x, T& xy_y) const
  138. {
  139. xy_x = asin(cos(lp_lat) * sin(lp_lon));
  140. xy_y = atan2(tan(lp_lat) , cos(lp_lon)) - this->m_par.phi0;
  141. }
  142. // INVERSE(s_inverse) spheroid
  143. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  144. inline void inv(T& xy_x, T& xy_y, T& lp_lon, T& lp_lat) const
  145. {
  146. T dd = xy_y + this->m_par.phi0;
  147. lp_lat = asin(sin(dd) * cos(xy_x));
  148. lp_lon = atan2(tan(xy_x), cos(dd));
  149. }
  150. static inline std::string get_name()
  151. {
  152. return "cass_spheroid";
  153. }
  154. };
  155. // Cassini
  156. template <typename Parameters, typename T>
  157. inline void setup_cass(Parameters& par, par_cass<T>& proj_parm)
  158. {
  159. if (par.es) {
  160. proj_parm.en = pj_enfn<T>(par.es);
  161. proj_parm.m0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en);
  162. } else {
  163. }
  164. }
  165. }} // namespace detail::cass
  166. #endif // doxygen
  167. /*!
  168. \brief Cassini projection
  169. \ingroup projections
  170. \tparam Geographic latlong point type
  171. \tparam Cartesian xy point type
  172. \tparam Parameters parameter type
  173. \par Projection characteristics
  174. - Cylindrical
  175. - Spheroid
  176. - Ellipsoid
  177. \par Example
  178. \image html ex_cass.gif
  179. */
  180. template <typename T, typename Parameters>
  181. struct cass_ellipsoid : public detail::cass::base_cass_ellipsoid<T, Parameters>
  182. {
  183. inline cass_ellipsoid(const Parameters& par) : detail::cass::base_cass_ellipsoid<T, Parameters>(par)
  184. {
  185. detail::cass::setup_cass(this->m_par, this->m_proj_parm);
  186. }
  187. };
  188. /*!
  189. \brief Cassini projection
  190. \ingroup projections
  191. \tparam Geographic latlong point type
  192. \tparam Cartesian xy point type
  193. \tparam Parameters parameter type
  194. \par Projection characteristics
  195. - Cylindrical
  196. - Spheroid
  197. - Ellipsoid
  198. \par Example
  199. \image html ex_cass.gif
  200. */
  201. template <typename T, typename Parameters>
  202. struct cass_spheroid : public detail::cass::base_cass_spheroid<T, Parameters>
  203. {
  204. inline cass_spheroid(const Parameters& par) : detail::cass::base_cass_spheroid<T, Parameters>(par)
  205. {
  206. detail::cass::setup_cass(this->m_par, this->m_proj_parm);
  207. }
  208. };
  209. #ifndef DOXYGEN_NO_DETAIL
  210. namespace detail
  211. {
  212. // Static projection
  213. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION(srs::par4::cass, cass_spheroid, cass_ellipsoid)
  214. // Factory entry(s)
  215. template <typename T, typename Parameters>
  216. class cass_entry : public detail::factory_entry<T, Parameters>
  217. {
  218. public :
  219. virtual base_v<T, Parameters>* create_new(const Parameters& par) const
  220. {
  221. if (par.es)
  222. return new base_v_fi<cass_ellipsoid<T, Parameters>, T, Parameters>(par);
  223. else
  224. return new base_v_fi<cass_spheroid<T, Parameters>, T, Parameters>(par);
  225. }
  226. };
  227. template <typename T, typename Parameters>
  228. inline void cass_init(detail::base_factory<T, Parameters>& factory)
  229. {
  230. factory.add_to_factory("cass", new cass_entry<T, Parameters>);
  231. }
  232. } // namespace detail
  233. #endif // doxygen
  234. } // namespace projections
  235. }} // namespace boost::geometry
  236. #endif // BOOST_GEOMETRY_PROJECTIONS_CASS_HPP