eqdc.hpp 8.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218
  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_EQDC_HPP
  31. #define BOOST_GEOMETRY_PROJECTIONS_EQDC_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_msfn.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. #include <boost/math/special_functions/hypot.hpp>
  41. namespace boost { namespace geometry
  42. {
  43. namespace projections
  44. {
  45. #ifndef DOXYGEN_NO_DETAIL
  46. namespace detail { namespace eqdc
  47. {
  48. static const double epsilon10 = 1.e-10;
  49. template <typename T>
  50. struct par_eqdc
  51. {
  52. T phi1;
  53. T phi2;
  54. T n;
  55. T rho0;
  56. T c;
  57. detail::en<T> en;
  58. bool ellips;
  59. };
  60. template <typename T, typename Parameters>
  61. struct base_eqdc_ellipsoid
  62. {
  63. par_eqdc<T> m_proj_parm;
  64. // FORWARD(e_forward) sphere & ellipsoid
  65. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  66. inline void fwd(Parameters const& , T lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
  67. {
  68. T rho = 0.0;
  69. rho = this->m_proj_parm.c - (this->m_proj_parm.ellips ? pj_mlfn(lp_lat, sin(lp_lat),
  70. cos(lp_lat), this->m_proj_parm.en) : lp_lat);
  71. xy_x = rho * sin( lp_lon *= this->m_proj_parm.n );
  72. xy_y = this->m_proj_parm.rho0 - rho * cos(lp_lon);
  73. }
  74. // INVERSE(e_inverse) sphere & ellipsoid
  75. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  76. inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
  77. {
  78. static T const half_pi = detail::half_pi<T>();
  79. T rho = 0.0;
  80. if ((rho = boost::math::hypot(xy_x, xy_y = this->m_proj_parm.rho0 - xy_y)) != 0.0 ) {
  81. if (this->m_proj_parm.n < 0.) {
  82. rho = -rho;
  83. xy_x = -xy_x;
  84. xy_y = -xy_y;
  85. }
  86. lp_lat = this->m_proj_parm.c - rho;
  87. if (this->m_proj_parm.ellips)
  88. lp_lat = pj_inv_mlfn(lp_lat, par.es, this->m_proj_parm.en);
  89. lp_lon = atan2(xy_x, xy_y) / this->m_proj_parm.n;
  90. } else {
  91. lp_lon = 0.;
  92. lp_lat = this->m_proj_parm.n > 0. ? half_pi : -half_pi;
  93. }
  94. }
  95. static inline std::string get_name()
  96. {
  97. return "eqdc_ellipsoid";
  98. }
  99. };
  100. // Equidistant Conic
  101. template <typename Params, typename Parameters, typename T>
  102. inline void setup_eqdc(Params const& params, Parameters& par, par_eqdc<T>& proj_parm)
  103. {
  104. T cosphi, sinphi;
  105. int secant;
  106. proj_parm.phi1 = pj_get_param_r<T, srs::spar::lat_1>(params, "lat_1", srs::dpar::lat_1);
  107. proj_parm.phi2 = pj_get_param_r<T, srs::spar::lat_2>(params, "lat_2", srs::dpar::lat_2);
  108. if (fabs(proj_parm.phi1 + proj_parm.phi2) < epsilon10)
  109. BOOST_THROW_EXCEPTION( projection_exception(error_conic_lat_equal) );
  110. proj_parm.en = pj_enfn<T>(par.es);
  111. proj_parm.n = sinphi = sin(proj_parm.phi1);
  112. cosphi = cos(proj_parm.phi1);
  113. secant = fabs(proj_parm.phi1 - proj_parm.phi2) >= epsilon10;
  114. if( (proj_parm.ellips = (par.es > 0.)) ) {
  115. double ml1, m1;
  116. m1 = pj_msfn(sinphi, cosphi, par.es);
  117. ml1 = pj_mlfn(proj_parm.phi1, sinphi, cosphi, proj_parm.en);
  118. if (secant) { /* secant cone */
  119. sinphi = sin(proj_parm.phi2);
  120. cosphi = cos(proj_parm.phi2);
  121. proj_parm.n = (m1 - pj_msfn(sinphi, cosphi, par.es)) /
  122. (pj_mlfn(proj_parm.phi2, sinphi, cosphi, proj_parm.en) - ml1);
  123. }
  124. proj_parm.c = ml1 + m1 / proj_parm.n;
  125. proj_parm.rho0 = proj_parm.c - pj_mlfn(par.phi0, sin(par.phi0),
  126. cos(par.phi0), proj_parm.en);
  127. } else {
  128. if (secant)
  129. proj_parm.n = (cosphi - cos(proj_parm.phi2)) / (proj_parm.phi2 - proj_parm.phi1);
  130. proj_parm.c = proj_parm.phi1 + cos(proj_parm.phi1) / proj_parm.n;
  131. proj_parm.rho0 = proj_parm.c - par.phi0;
  132. }
  133. }
  134. }} // namespace detail::eqdc
  135. #endif // doxygen
  136. /*!
  137. \brief Equidistant Conic projection
  138. \ingroup projections
  139. \tparam Geographic latlong point type
  140. \tparam Cartesian xy point type
  141. \tparam Parameters parameter type
  142. \par Projection characteristics
  143. - Conic
  144. - Spheroid
  145. - Ellipsoid
  146. \par Projection parameters
  147. - lat_1: Latitude of first standard parallel (degrees)
  148. - lat_2: Latitude of second standard parallel (degrees)
  149. \par Example
  150. \image html ex_eqdc.gif
  151. */
  152. template <typename T, typename Parameters>
  153. struct eqdc_ellipsoid : public detail::eqdc::base_eqdc_ellipsoid<T, Parameters>
  154. {
  155. template <typename Params>
  156. inline eqdc_ellipsoid(Params const& params, Parameters const& par)
  157. {
  158. detail::eqdc::setup_eqdc(params, par, this->m_proj_parm);
  159. }
  160. };
  161. #ifndef DOXYGEN_NO_DETAIL
  162. namespace detail
  163. {
  164. // Static projection
  165. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_eqdc, eqdc_ellipsoid)
  166. // Factory entry(s)
  167. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(eqdc_entry, eqdc_ellipsoid)
  168. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(eqdc_init)
  169. {
  170. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(eqdc, eqdc_entry);
  171. }
  172. } // namespace detail
  173. #endif // doxygen
  174. } // namespace projections
  175. }} // namespace boost::geometry
  176. #endif // BOOST_GEOMETRY_PROJECTIONS_EQDC_HPP