bonne.hpp 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265
  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_BONNE_HPP
  31. #define BOOST_GEOMETRY_PROJECTIONS_BONNE_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/factory_entry.hpp>
  35. #include <boost/geometry/srs/projections/impl/pj_mlfn.hpp>
  36. #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  37. #include <boost/geometry/srs/projections/impl/projects.hpp>
  38. #include <boost/geometry/util/math.hpp>
  39. #include <boost/math/special_functions/hypot.hpp>
  40. namespace boost { namespace geometry
  41. {
  42. namespace projections
  43. {
  44. #ifndef DOXYGEN_NO_DETAIL
  45. namespace detail { namespace bonne
  46. {
  47. static const double epsilon10 = 1e-10;
  48. template <typename T>
  49. struct par_bonne
  50. {
  51. T phi1;
  52. T cphi1;
  53. T am1;
  54. T m1;
  55. detail::en<T> en;
  56. };
  57. template <typename T, typename Parameters>
  58. struct base_bonne_ellipsoid
  59. {
  60. par_bonne<T> m_proj_parm;
  61. // FORWARD(e_forward) ellipsoid
  62. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  63. inline void fwd(Parameters const& par, T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  64. {
  65. T rh, E, c;
  66. rh = this->m_proj_parm.am1 + this->m_proj_parm.m1 - pj_mlfn(lp_lat, E = sin(lp_lat), c = cos(lp_lat), this->m_proj_parm.en);
  67. E = c * lp_lon / (rh * sqrt(1. - par.es * E * E));
  68. xy_x = rh * sin(E);
  69. xy_y = this->m_proj_parm.am1 - rh * cos(E);
  70. }
  71. // INVERSE(e_inverse) ellipsoid
  72. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  73. inline void inv(Parameters const& par, T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
  74. {
  75. static const T half_pi = detail::half_pi<T>();
  76. T s, rh;
  77. rh = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.am1 - xy_y);
  78. lp_lat = pj_inv_mlfn(this->m_proj_parm.am1 + this->m_proj_parm.m1 - rh, par.es, this->m_proj_parm.en);
  79. if ((s = fabs(lp_lat)) < half_pi) {
  80. s = sin(lp_lat);
  81. lp_lon = rh * atan2(xy_x, xy_y) *
  82. sqrt(1. - par.es * s * s) / cos(lp_lat);
  83. } else if (fabs(s - half_pi) <= epsilon10)
  84. lp_lon = 0.;
  85. else
  86. BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
  87. }
  88. static inline std::string get_name()
  89. {
  90. return "bonne_ellipsoid";
  91. }
  92. };
  93. template <typename T, typename Parameters>
  94. struct base_bonne_spheroid
  95. {
  96. par_bonne<T> m_proj_parm;
  97. // FORWARD(s_forward) spheroid
  98. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  99. inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  100. {
  101. T E, rh;
  102. rh = this->m_proj_parm.cphi1 + this->m_proj_parm.phi1 - lp_lat;
  103. if (fabs(rh) > epsilon10) {
  104. xy_x = rh * sin(E = lp_lon * cos(lp_lat) / rh);
  105. xy_y = this->m_proj_parm.cphi1 - rh * cos(E);
  106. } else
  107. xy_x = xy_y = 0.;
  108. }
  109. // INVERSE(s_inverse) spheroid
  110. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  111. inline void inv(Parameters const& , T const& xy_x, T xy_y, T& lp_lon, T& lp_lat) const
  112. {
  113. static const T half_pi = detail::half_pi<T>();
  114. T rh;
  115. rh = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.cphi1 - xy_y);
  116. lp_lat = this->m_proj_parm.cphi1 + this->m_proj_parm.phi1 - rh;
  117. if (fabs(lp_lat) > half_pi) {
  118. BOOST_THROW_EXCEPTION( projection_exception(error_tolerance_condition) );
  119. }
  120. if (fabs(fabs(lp_lat) - half_pi) <= epsilon10)
  121. lp_lon = 0.;
  122. else
  123. lp_lon = rh * atan2(xy_x, xy_y) / cos(lp_lat);
  124. }
  125. static inline std::string get_name()
  126. {
  127. return "bonne_spheroid";
  128. }
  129. };
  130. // Bonne (Werner lat_1=90)
  131. template <typename Params, typename Parameters, typename T>
  132. inline void setup_bonne(Params const& params, Parameters const& par, par_bonne<T>& proj_parm)
  133. {
  134. static const T half_pi = detail::half_pi<T>();
  135. T c;
  136. proj_parm.phi1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
  137. if (fabs(proj_parm.phi1) < epsilon10)
  138. BOOST_THROW_EXCEPTION( projection_exception(error_lat1_is_zero) );
  139. if (par.es != 0.0) {
  140. proj_parm.en = pj_enfn<T>(par.es);
  141. proj_parm.m1 = pj_mlfn(proj_parm.phi1, proj_parm.am1 = sin(proj_parm.phi1),
  142. c = cos(proj_parm.phi1), proj_parm.en);
  143. proj_parm.am1 = c / (sqrt(1. - par.es * proj_parm.am1 * proj_parm.am1) * proj_parm.am1);
  144. } else {
  145. if (fabs(proj_parm.phi1) + epsilon10 >= half_pi)
  146. proj_parm.cphi1 = 0.;
  147. else
  148. proj_parm.cphi1 = 1. / tan(proj_parm.phi1);
  149. }
  150. }
  151. }} // namespace detail::bonne
  152. #endif // doxygen
  153. /*!
  154. \brief Bonne (Werner lat_1=90) 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. - Conic
  161. - Spheroid
  162. - Ellipsoid
  163. \par Projection parameters
  164. - lat_1: Latitude of first standard parallel (degrees)
  165. \par Example
  166. \image html ex_bonne.gif
  167. */
  168. template <typename T, typename Parameters>
  169. struct bonne_ellipsoid : public detail::bonne::base_bonne_ellipsoid<T, Parameters>
  170. {
  171. template <typename Params>
  172. inline bonne_ellipsoid(Params const& params, Parameters const& par)
  173. {
  174. detail::bonne::setup_bonne(params, par, this->m_proj_parm);
  175. }
  176. };
  177. /*!
  178. \brief Bonne (Werner lat_1=90) projection
  179. \ingroup projections
  180. \tparam Geographic latlong point type
  181. \tparam Cartesian xy point type
  182. \tparam Parameters parameter type
  183. \par Projection characteristics
  184. - Conic
  185. - Spheroid
  186. - Ellipsoid
  187. \par Projection parameters
  188. - lat_1: Latitude of first standard parallel (degrees)
  189. \par Example
  190. \image html ex_bonne.gif
  191. */
  192. template <typename T, typename Parameters>
  193. struct bonne_spheroid : public detail::bonne::base_bonne_spheroid<T, Parameters>
  194. {
  195. template <typename Params>
  196. inline bonne_spheroid(Params const& params, Parameters const& par)
  197. {
  198. detail::bonne::setup_bonne(params, par, this->m_proj_parm);
  199. }
  200. };
  201. #ifndef DOXYGEN_NO_DETAIL
  202. namespace detail
  203. {
  204. // Static projection
  205. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_bonne, bonne_spheroid, bonne_ellipsoid)
  206. // Factory entry(s)
  207. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(bonne_entry, bonne_spheroid, bonne_ellipsoid)
  208. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(bonne_init)
  209. {
  210. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(bonne, bonne_entry);
  211. }
  212. } // namespace detail
  213. #endif // doxygen
  214. } // namespace projections
  215. }} // namespace boost::geometry
  216. #endif // BOOST_GEOMETRY_PROJECTIONS_BONNE_HPP