cass.hpp 10 KB

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