nzmg.hpp 8.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230
  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. // Purpose: Implementation of the nzmg (New Zealand Map Grid) projection.
  16. // Very loosely based upon DMA code by Bradford W. Drew
  17. // Author: Gerald Evenden
  18. // Copyright (c) 1995, Gerald Evenden
  19. // Permission is hereby granted, free of charge, to any person obtaining a
  20. // copy of this software and associated documentation files (the "Software"),
  21. // to deal in the Software without restriction, including without limitation
  22. // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  23. // and/or sell copies of the Software, and to permit persons to whom the
  24. // Software is furnished to do so, subject to the following conditions:
  25. // The above copyright notice and this permission notice shall be included
  26. // in all copies or substantial portions of the Software.
  27. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  28. // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  29. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  30. // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  31. // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  32. // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  33. // DEALINGS IN THE SOFTWARE.
  34. #ifndef BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
  35. #define BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
  36. #include <boost/geometry/util/math.hpp>
  37. #include <boost/geometry/srs/projections/impl/base_static.hpp>
  38. #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  39. #include <boost/geometry/srs/projections/impl/projects.hpp>
  40. #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  41. #include <boost/geometry/srs/projections/impl/pj_zpoly1.hpp>
  42. namespace boost { namespace geometry
  43. {
  44. namespace projections
  45. {
  46. #ifndef DOXYGEN_NO_DETAIL
  47. namespace detail { namespace nzmg
  48. {
  49. static const double epsilon = 1e-10;
  50. static const int Nbf = 5;
  51. static const int Ntpsi = 9;
  52. static const int Ntphi = 8;
  53. template <typename T>
  54. inline T sec5_to_rad() { return 0.4848136811095359935899141023; }
  55. template <typename T>
  56. inline T rad_to_sec5() { return 2.062648062470963551564733573; }
  57. template <typename T>
  58. inline const pj_complex<T> * bf()
  59. {
  60. static const pj_complex<T> result[] = {
  61. {.7557853228, 0.0},
  62. {.249204646, .003371507},
  63. {-.001541739, .041058560},
  64. {-.10162907, .01727609},
  65. {-.26623489, -.36249218},
  66. {-.6870983, -1.1651967}
  67. };
  68. return result;
  69. }
  70. template <typename T>
  71. inline const T * tphi()
  72. {
  73. static const T result[] = { 1.5627014243, .5185406398, -.03333098,
  74. -.1052906, -.0368594, .007317,
  75. .01220, .00394, -.0013 };
  76. return result;
  77. }
  78. template <typename T>
  79. inline const T * tpsi()
  80. {
  81. static const T result[] = { .6399175073, -.1358797613, .063294409, -.02526853, .0117879,
  82. -.0055161, .0026906, -.001333, .00067, -.00034 };
  83. return result;
  84. }
  85. template <typename T, typename Parameters>
  86. struct base_nzmg_ellipsoid
  87. {
  88. // FORWARD(e_forward) ellipsoid
  89. // Project coordinates from geographic (lon, lat) to cartesian (x, y)
  90. inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
  91. {
  92. static const T rad_to_sec5 = nzmg::rad_to_sec5<T>();
  93. pj_complex<T> p;
  94. const T * C;
  95. int i;
  96. lp_lat = (lp_lat - par.phi0) * rad_to_sec5;
  97. for (p.r = *(C = tpsi<T>() + (i = Ntpsi)); i ; --i)
  98. p.r = *--C + lp_lat * p.r;
  99. p.r *= lp_lat;
  100. p.i = lp_lon;
  101. p = pj_zpoly1(p, bf<T>(), Nbf);
  102. xy_x = p.i;
  103. xy_y = p.r;
  104. }
  105. // INVERSE(e_inverse) ellipsoid
  106. // Project coordinates from cartesian (x, y) to geographic (lon, lat)
  107. inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
  108. {
  109. static const T sec5_to_rad = nzmg::sec5_to_rad<T>();
  110. int nn, i;
  111. pj_complex<T> p, f, fp, dp;
  112. T den;
  113. const T* C;
  114. p.r = xy_y;
  115. p.i = xy_x;
  116. for (nn = 20; nn ;--nn) {
  117. f = pj_zpolyd1(p, bf<T>(), Nbf, &fp);
  118. f.r -= xy_y;
  119. f.i -= xy_x;
  120. den = fp.r * fp.r + fp.i * fp.i;
  121. p.r += dp.r = -(f.r * fp.r + f.i * fp.i) / den;
  122. p.i += dp.i = -(f.i * fp.r - f.r * fp.i) / den;
  123. if ((fabs(dp.r) + fabs(dp.i)) <= epsilon)
  124. break;
  125. }
  126. if (nn) {
  127. lp_lon = p.i;
  128. for (lp_lat = *(C = tphi<T>() + (i = Ntphi)); i ; --i)
  129. lp_lat = *--C + p.r * lp_lat;
  130. lp_lat = par.phi0 + p.r * lp_lat * sec5_to_rad;
  131. } else
  132. lp_lon = lp_lat = HUGE_VAL;
  133. }
  134. static inline std::string get_name()
  135. {
  136. return "nzmg_ellipsoid";
  137. }
  138. };
  139. // New Zealand Map Grid
  140. template <typename Parameters>
  141. inline void setup_nzmg(Parameters& par)
  142. {
  143. typedef typename Parameters::type calc_t;
  144. static const calc_t d2r = geometry::math::d2r<calc_t>();
  145. /* force to International major axis */
  146. par.a = 6378388.0;
  147. par.ra = 1. / par.a;
  148. par.lam0 = 173. * d2r;
  149. par.phi0 = -41. * d2r;
  150. par.x0 = 2510000.;
  151. par.y0 = 6023150.;
  152. }
  153. }} // namespace detail::nzmg
  154. #endif // doxygen
  155. /*!
  156. \brief New Zealand Map Grid projection
  157. \ingroup projections
  158. \tparam Geographic latlong point type
  159. \tparam Cartesian xy point type
  160. \tparam Parameters parameter type
  161. \par Projection characteristics
  162. - Fixed Earth
  163. \par Example
  164. \image html ex_nzmg.gif
  165. */
  166. template <typename T, typename Parameters>
  167. struct nzmg_ellipsoid : public detail::nzmg::base_nzmg_ellipsoid<T, Parameters>
  168. {
  169. template <typename Params>
  170. inline nzmg_ellipsoid(Params const& , Parameters & par)
  171. {
  172. detail::nzmg::setup_nzmg(par);
  173. }
  174. };
  175. #ifndef DOXYGEN_NO_DETAIL
  176. namespace detail
  177. {
  178. // Static projection
  179. BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_nzmg, nzmg_ellipsoid)
  180. // Factory entry(s)
  181. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(nzmg_entry, nzmg_ellipsoid)
  182. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(nzmg_init)
  183. {
  184. BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(nzmg, nzmg_entry)
  185. }
  186. } // namespace detail
  187. #endif // doxygen
  188. } // namespace projections
  189. }} // namespace boost::geometry
  190. #endif // BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP