cass.hpp 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255
  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_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 <typename T, typename Parameters>
  67. struct base_cass_ellipsoid
  68. {
  69. par_cass<T> m_proj_parm;
  70. // FORWARD(e_forward) ellipsoid
  71. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  72. inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  73. {
  74. static const T C1 = cass::C1<T>();
  75. static const T C2 = cass::C2<T>();
  76. static const T C3 = cass::C3<T>();
  77. T n = sin(lp_lat);
  78. T c = cos(lp_lat);
  79. xy_y = pj_mlfn(lp_lat, n, c, this->m_proj_parm.en);
  80. n = 1./sqrt(1. - par.es * n * n);
  81. T tn = tan(lp_lat);
  82. T t = tn * tn;
  83. T a1 = lp_lon * c;
  84. c *= par.es * c / (1 - par.es);
  85. T a2 = a1 * a1;
  86. xy_x = n * a1 * (1. - a2 * t *
  87. (C1 - (8. - t + 8. * c) * a2 * C2));
  88. xy_y -= this->m_proj_parm.m0 - n * tn * a2 *
  89. (.5 + (5. - t + 6. * c) * a2 * C3);
  90. }
  91. // INVERSE(e_inverse) ellipsoid
  92. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  93. inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
  94. {
  95. static const T C3 = cass::C3<T>();
  96. static const T C4 = cass::C4<T>();
  97. static const T C5 = cass::C5<T>();
  98. T ph1;
  99. ph1 = pj_inv_mlfn(this->m_proj_parm.m0 + xy_y, par.es, this->m_proj_parm.en);
  100. T tn = tan(ph1); T t = tn * tn;
  101. T n = sin(ph1);
  102. T r = 1. / (1. - par.es * n * n);
  103. n = sqrt(r);
  104. r *= (1. - par.es) * n;
  105. T dd = xy_x / n;
  106. T d2 = dd * dd;
  107. lp_lat = ph1 - (n * tn / r) * d2 *
  108. (.5 - (1. + 3. * t) * d2 * C3);
  109. lp_lon = dd * (1. + t * d2 *
  110. (-C4 + (1. + 3. * t) * d2 * C5)) / cos(ph1);
  111. }
  112. static inline std::string get_name()
  113. {
  114. return "cass_ellipsoid";
  115. }
  116. };
  117. template <typename T, typename Parameters>
  118. struct base_cass_spheroid
  119. {
  120. par_cass<T> m_proj_parm;
  121. // FORWARD(s_forward) spheroid
  122. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  123. inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  124. {
  125. xy_x = asin(cos(lp_lat) * sin(lp_lon));
  126. xy_y = atan2(tan(lp_lat) , cos(lp_lon)) - par.phi0;
  127. }
  128. // INVERSE(s_inverse) spheroid
  129. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  130. inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
  131. {
  132. T dd = xy_y + par.phi0;
  133. lp_lat = asin(sin(dd) * cos(xy_x));
  134. lp_lon = atan2(tan(xy_x), cos(dd));
  135. }
  136. static inline std::string get_name()
  137. {
  138. return "cass_spheroid";
  139. }
  140. };
  141. // Cassini
  142. template <typename Parameters, typename T>
  143. inline void setup_cass(Parameters& par, par_cass<T>& proj_parm)
  144. {
  145. if (par.es) {
  146. proj_parm.en = pj_enfn<T>(par.es);
  147. proj_parm.m0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en);
  148. } else {
  149. }
  150. }
  151. }} // namespace detail::cass
  152. #endif // doxygen
  153. /*!
  154. \brief Cassini projection
  155. \ingroup projections
  156. \tparam Geographic latlong point type
  157. \tparam Cartesian xy point type
  158. \tparam Parameters parameter type
  159. \par Projection characteristics
  160. - Cylindrical
  161. - Spheroid
  162. - Ellipsoid
  163. \par Example
  164. \image html ex_cass.gif
  165. */
  166. template <typename T, typename Parameters>
  167. struct cass_ellipsoid : public detail::cass::base_cass_ellipsoid<T, Parameters>
  168. {
  169. template <typename Params>
  170. inline cass_ellipsoid(Params const& , Parameters const& par)
  171. {
  172. detail::cass::setup_cass(par, this->m_proj_parm);
  173. }
  174. };
  175. /*!
  176. \brief Cassini projection
  177. \ingroup projections
  178. \tparam Geographic latlong point type
  179. \tparam Cartesian xy point type
  180. \tparam Parameters parameter type
  181. \par Projection characteristics
  182. - Cylindrical
  183. - Spheroid
  184. - Ellipsoid
  185. \par Example
  186. \image html ex_cass.gif
  187. */
  188. template <typename T, typename Parameters>
  189. struct cass_spheroid : public detail::cass::base_cass_spheroid<T, Parameters>
  190. {
  191. template <typename Params>
  192. inline cass_spheroid(Params const& , Parameters const& par)
  193. {
  194. detail::cass::setup_cass(par, this->m_proj_parm);
  195. }
  196. };
  197. #ifndef DOXYGEN_NO_DETAIL
  198. namespace detail
  199. {
  200. // Static projection
  201. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_cass, cass_spheroid, cass_ellipsoid)
  202. // Factory entry(s)
  203. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(cass_entry, cass_spheroid, cass_ellipsoid)
  204. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(cass_init)
  205. {
  206. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(cass, cass_entry);
  207. }
  208. } // namespace detail
  209. #endif // doxygen
  210. } // namespace projections
  211. }} // namespace boost::geometry
  212. #endif // BOOST_GEOMETRY_PROJECTIONS_CASS_HPP