TransverseMercator.hpp 9.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206
  1. /**
  2. * \file TransverseMercator.hpp
  3. * \brief Header for GeographicLib::TransverseMercator class
  4. *
  5. * Copyright (c) Charles Karney (2008-2020) <charles@karney.com> and licensed
  6. * under the MIT/X11 License. For more information, see
  7. * https://geographiclib.sourceforge.io/
  8. **********************************************************************/
  9. #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP)
  10. #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP 1
  11. #include <GeographicLib/Constants.hpp>
  12. #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER)
  13. /**
  14. * The order of the series approximation used in TransverseMercator.
  15. * GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER can be set to any integer in [4, 8].
  16. **********************************************************************/
  17. # define GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER \
  18. (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \
  19. (GEOGRAPHICLIB_PRECISION == 1 ? 4 : 8))
  20. #endif
  21. namespace GeographicLib {
  22. /**
  23. * \brief Transverse Mercator projection
  24. *
  25. * This uses Kr&uuml;ger's method which evaluates the projection and its
  26. * inverse in terms of a series. See
  27. * - L. Kr&uuml;ger,
  28. * <a href="https://doi.org/10.2312/GFZ.b103-krueger28"> Konforme
  29. * Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the
  30. * ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New
  31. * Series 52, 172 pp. (1912).
  32. * - C. F. F. Karney,
  33. * <a href="https://doi.org/10.1007/s00190-011-0445-3">
  34. * Transverse Mercator with an accuracy of a few nanometers,</a>
  35. * J. Geodesy 85(8), 475--485 (Aug. 2011);
  36. * preprint
  37. * <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>.
  38. *
  39. * Kr&uuml;ger's method has been extended from 4th to 6th order. The maximum
  40. * error is 5 nm (5 nanometers), ground distance, for all positions within 35
  41. * degrees of the central meridian. The error in the convergence is 2
  42. * &times; 10<sup>&minus;15</sup>&quot; and the relative error in the scale
  43. * is 6 &times; 10<sup>&minus;12</sup>%%. See Sec. 4 of
  44. * <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for details.
  45. * The speed penalty in going to 6th order is only about 1%.
  46. *
  47. * There's a singularity in the projection at &phi; = 0&deg;, &lambda;
  48. * &minus; &lambda;<sub>0</sub> = &plusmn;(1 &minus; \e e)90&deg; (&asymp;
  49. * &plusmn;82.6&deg; for the WGS84 ellipsoid), where \e e is the
  50. * eccentricity. Beyond this point, the series ceases to converge and the
  51. * results from this method will be garbage. To be on the safe side, don't
  52. * use this method if the angular distance from the central meridian exceeds
  53. * (1 &minus; 2e)90&deg; (&asymp; 75&deg; for the WGS84 ellipsoid)
  54. *
  55. * TransverseMercatorExact is an alternative implementation of the projection
  56. * using exact formulas which yield accurate (to 8 nm) results over the
  57. * entire ellipsoid.
  58. *
  59. * The ellipsoid parameters and the central scale are set in the constructor.
  60. * The central meridian (which is a trivial shift of the longitude) is
  61. * specified as the \e lon0 argument of the TransverseMercator::Forward and
  62. * TransverseMercator::Reverse functions. The latitude of origin is taken to
  63. * be the equator. There is no provision in this class for specifying a
  64. * false easting or false northing or a different latitude of origin.
  65. * However these are can be simply included by the calling function. For
  66. * example, the UTMUPS class applies the false easting and false northing for
  67. * the UTM projections. A more complicated example is the British National
  68. * Grid (<a href="https://www.spatialreference.org/ref/epsg/7405/">
  69. * EPSG:7405</a>) which requires the use of a latitude of origin. This is
  70. * implemented by the GeographicLib::OSGB class.
  71. *
  72. * This class also returns the meridian convergence \e gamma and scale \e k.
  73. * The meridian convergence is the bearing of grid north (the \e y axis)
  74. * measured clockwise from true north.
  75. *
  76. * See TransverseMercator.cpp for more information on the implementation.
  77. *
  78. * See \ref transversemercator for a discussion of this projection.
  79. *
  80. * Example of use:
  81. * \include example-TransverseMercator.cpp
  82. *
  83. * <a href="TransverseMercatorProj.1.html">TransverseMercatorProj</a> is a
  84. * command-line utility providing access to the functionality of
  85. * TransverseMercator and TransverseMercatorExact.
  86. **********************************************************************/
  87. class GEOGRAPHICLIB_EXPORT TransverseMercator {
  88. private:
  89. typedef Math::real real;
  90. static const int maxpow_ = GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER;
  91. static const int numit_ = 5;
  92. real _a, _f, _k0, _e2, _es, _e2m, _c, _n;
  93. // _alp[0] and _bet[0] unused
  94. real _a1, _b1, _alp[maxpow_ + 1], _bet[maxpow_ + 1];
  95. friend class Ellipsoid; // For access to taupf, tauf.
  96. public:
  97. /**
  98. * Constructor for a ellipsoid with
  99. *
  100. * @param[in] a equatorial radius (meters).
  101. * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
  102. * Negative \e f gives a prolate ellipsoid.
  103. * @param[in] k0 central scale factor.
  104. * @exception GeographicErr if \e a, (1 &minus; \e f) \e a, or \e k0 is
  105. * not positive.
  106. **********************************************************************/
  107. TransverseMercator(real a, real f, real k0);
  108. /**
  109. * Forward projection, from geographic to transverse Mercator.
  110. *
  111. * @param[in] lon0 central meridian of the projection (degrees).
  112. * @param[in] lat latitude of point (degrees).
  113. * @param[in] lon longitude of point (degrees).
  114. * @param[out] x easting of point (meters).
  115. * @param[out] y northing of point (meters).
  116. * @param[out] gamma meridian convergence at point (degrees).
  117. * @param[out] k scale of projection at point.
  118. *
  119. * No false easting or northing is added. \e lat should be in the range
  120. * [&minus;90&deg;, 90&deg;].
  121. **********************************************************************/
  122. void Forward(real lon0, real lat, real lon,
  123. real& x, real& y, real& gamma, real& k) const;
  124. /**
  125. * Reverse projection, from transverse Mercator to geographic.
  126. *
  127. * @param[in] lon0 central meridian of the projection (degrees).
  128. * @param[in] x easting of point (meters).
  129. * @param[in] y northing of point (meters).
  130. * @param[out] lat latitude of point (degrees).
  131. * @param[out] lon longitude of point (degrees).
  132. * @param[out] gamma meridian convergence at point (degrees).
  133. * @param[out] k scale of projection at point.
  134. *
  135. * No false easting or northing is added. The value of \e lon returned is
  136. * in the range [&minus;180&deg;, 180&deg;].
  137. **********************************************************************/
  138. void Reverse(real lon0, real x, real y,
  139. real& lat, real& lon, real& gamma, real& k) const;
  140. /**
  141. * TransverseMercator::Forward without returning the convergence and scale.
  142. **********************************************************************/
  143. void Forward(real lon0, real lat, real lon,
  144. real& x, real& y) const {
  145. real gamma, k;
  146. Forward(lon0, lat, lon, x, y, gamma, k);
  147. }
  148. /**
  149. * TransverseMercator::Reverse without returning the convergence and scale.
  150. **********************************************************************/
  151. void Reverse(real lon0, real x, real y,
  152. real& lat, real& lon) const {
  153. real gamma, k;
  154. Reverse(lon0, x, y, lat, lon, gamma, k);
  155. }
  156. /** \name Inspector functions
  157. **********************************************************************/
  158. ///@{
  159. /**
  160. * @return \e a the equatorial radius of the ellipsoid (meters). This is
  161. * the value used in the constructor.
  162. **********************************************************************/
  163. Math::real EquatorialRadius() const { return _a; }
  164. /**
  165. * @return \e f the flattening of the ellipsoid. This is the value used in
  166. * the constructor.
  167. **********************************************************************/
  168. Math::real Flattening() const { return _f; }
  169. /**
  170. * @return \e k0 central scale for the projection. This is the value of \e
  171. * k0 used in the constructor and is the scale on the central meridian.
  172. **********************************************************************/
  173. Math::real CentralScale() const { return _k0; }
  174. /**
  175. * \deprecated An old name for EquatorialRadius().
  176. **********************************************************************/
  177. GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
  178. Math::real MajorRadius() const { return EquatorialRadius(); }
  179. ///@}
  180. /**
  181. * A global instantiation of TransverseMercator with the WGS84 ellipsoid
  182. * and the UTM scale factor. However, unlike UTM, no false easting or
  183. * northing is added.
  184. **********************************************************************/
  185. static const TransverseMercator& UTM();
  186. };
  187. } // namespace GeographicLib
  188. #endif // GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP