00001 /** 00002 * \file TransverseMercator.hpp 00003 * \brief Header for GeographicLib::TransverseMercator class 00004 * 00005 * Copyright (c) Charles Karney (2008-2014) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP) 00011 #define GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 00015 #if !defined(GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER) 00016 /** 00017 * The order of the series approximation used in TransverseMercator. 00018 * GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER can be set to any integer in [4, 8]. 00019 **********************************************************************/ 00020 # define GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER \ 00021 (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \ 00022 (GEOGRAPHICLIB_PRECISION == 1 ? 4 : 8)) 00023 #endif 00024 00025 namespace GeographicLib { 00026 00027 /** 00028 * \brief Transverse Mercator projection 00029 * 00030 * This uses Krüger's method which evaluates the projection and its 00031 * inverse in terms of a series. See 00032 * - L. Krüger, 00033 * <a href="http://dx.doi.org/10.2312/GFZ.b103-krueger28"> Konforme 00034 * Abbildung des Erdellipsoids in der Ebene</a> (Conformal mapping of the 00035 * ellipsoidal earth to the plane), Royal Prussian Geodetic Institute, New 00036 * Series 52, 172 pp. (1912). 00037 * - C. F. F. Karney, 00038 * <a href="http://dx.doi.org/10.1007/s00190-011-0445-3"> 00039 * Transverse Mercator with an accuracy of a few nanometers,</a> 00040 * J. Geodesy 85(8), 475--485 (Aug. 2011); 00041 * preprint 00042 * <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>. 00043 * 00044 * Krüger's method has been extended from 4th to 6th order. The maximum 00045 * error is 5 nm (5 nanometers), ground distance, for all positions within 35 00046 * degrees of the central meridian. The error in the convergence is 2 00047 * × 10<sup>−15</sup>" and the relative error in the scale 00048 * is 6 − 10<sup>−12</sup>%%. See Sec. 4 of 00049 * <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for details. 00050 * The speed penalty in going to 6th order is only about 1%. 00051 * TransverseMercatorExact is an alternative implementation of the projection 00052 * using exact formulas which yield accurate (to 8 nm) results over the 00053 * entire ellipsoid. 00054 * 00055 * The ellipsoid parameters and the central scale are set in the constructor. 00056 * The central meridian (which is a trivial shift of the longitude) is 00057 * specified as the \e lon0 argument of the TransverseMercator::Forward and 00058 * TransverseMercator::Reverse functions. The latitude of origin is taken to 00059 * be the equator. There is no provision in this class for specifying a 00060 * false easting or false northing or a different latitude of origin. 00061 * However these are can be simply included by the calling function. For 00062 * example, the UTMUPS class applies the false easting and false northing for 00063 * the UTM projections. A more complicated example is the British National 00064 * Grid (<a href="http://www.spatialreference.org/ref/epsg/7405/"> 00065 * EPSG:7405</a>) which requires the use of a latitude of origin. This is 00066 * implemented by the GeographicLib::OSGB class. 00067 * 00068 * See TransverseMercator.cpp for more information on the implementation. 00069 * 00070 * See \ref transversemercator for a discussion of this projection. 00071 * 00072 * Example of use: 00073 * \include example-TransverseMercator.cpp 00074 * 00075 * <a href="TransverseMercatorProj.1.html">TransverseMercatorProj</a> is a 00076 * command-line utility providing access to the functionality of 00077 * TransverseMercator and TransverseMercatorExact. 00078 **********************************************************************/ 00079 00080 class GEOGRAPHICLIB_EXPORT TransverseMercator { 00081 private: 00082 typedef Math::real real; 00083 static const int maxpow_ = GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER; 00084 static const int numit_ = 5; 00085 real tol_; 00086 real _a, _f, _k0, _e2, _e, _e2m, _c, _n; 00087 // _alp[0] and _bet[0] unused 00088 real _a1, _b1, _alp[maxpow_ + 1], _bet[maxpow_ + 1]; 00089 static inline real overflow() { 00090 // Overflow value s.t. atan(overflow_) = pi/2 00091 static const real 00092 overflow = 1 / Math::sq(std::numeric_limits<real>::epsilon()); 00093 return overflow; 00094 } 00095 // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right 00096 static inline real tanx(real x) { 00097 using std::tan; 00098 real t = tan(x); 00099 // Write the tests this way to ensure that tanx(NaN()) is NaN() 00100 return x >= 0 ? 00101 (!(t < 0) ? t : overflow()) : 00102 (!(t >= 0) ? t : -overflow()); 00103 } 00104 // Return e * atanh(e * x) for f >= 0, else return 00105 // - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 00106 inline real eatanhe(real x) const { 00107 using std::atan; 00108 return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * atan(_e * x); 00109 } 00110 real taupf(real tau) const; 00111 real tauf(real taup) const; 00112 00113 friend class Ellipsoid; // For access to taupf, tauf. 00114 public: 00115 00116 /** 00117 * Constructor for a ellipsoid with 00118 * 00119 * @param[in] a equatorial radius (meters). 00120 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00121 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set 00122 * flattening to 1/\e f. 00123 * @param[in] k0 central scale factor. 00124 * @exception GeographicErr if \e a, (1 − \e f) \e a, or \e k0 is 00125 * not positive. 00126 **********************************************************************/ 00127 TransverseMercator(real a, real f, real k0); 00128 00129 /** 00130 * Forward projection, from geographic to transverse Mercator. 00131 * 00132 * @param[in] lon0 central meridian of the projection (degrees). 00133 * @param[in] lat latitude of point (degrees). 00134 * @param[in] lon longitude of point (degrees). 00135 * @param[out] x easting of point (meters). 00136 * @param[out] y northing of point (meters). 00137 * @param[out] gamma meridian convergence at point (degrees). 00138 * @param[out] k scale of projection at point. 00139 * 00140 * No false easting or northing is added. \e lat should be in the range 00141 * [−90°, 90°]; \e lon and \e lon0 should be in the 00142 * range [−540°, 540°). 00143 **********************************************************************/ 00144 void Forward(real lon0, real lat, real lon, 00145 real& x, real& y, real& gamma, real& k) const; 00146 00147 /** 00148 * Reverse projection, from transverse Mercator to geographic. 00149 * 00150 * @param[in] lon0 central meridian of the projection (degrees). 00151 * @param[in] x easting of point (meters). 00152 * @param[in] y northing of point (meters). 00153 * @param[out] lat latitude of point (degrees). 00154 * @param[out] lon longitude of point (degrees). 00155 * @param[out] gamma meridian convergence at point (degrees). 00156 * @param[out] k scale of projection at point. 00157 * 00158 * No false easting or northing is added. \e lon0 should be in the range 00159 * [−540°, 540°). The value of \e lon returned is in 00160 * the range [−180°, 180°). 00161 **********************************************************************/ 00162 void Reverse(real lon0, real x, real y, 00163 real& lat, real& lon, real& gamma, real& k) const; 00164 00165 /** 00166 * TransverseMercator::Forward without returning the convergence and scale. 00167 **********************************************************************/ 00168 void Forward(real lon0, real lat, real lon, 00169 real& x, real& y) const { 00170 real gamma, k; 00171 Forward(lon0, lat, lon, x, y, gamma, k); 00172 } 00173 00174 /** 00175 * TransverseMercator::Reverse without returning the convergence and scale. 00176 **********************************************************************/ 00177 void Reverse(real lon0, real x, real y, 00178 real& lat, real& lon) const { 00179 real gamma, k; 00180 Reverse(lon0, x, y, lat, lon, gamma, k); 00181 } 00182 00183 /** \name Inspector functions 00184 **********************************************************************/ 00185 ///@{ 00186 /** 00187 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00188 * the value used in the constructor. 00189 **********************************************************************/ 00190 Math::real MajorRadius() const { return _a; } 00191 00192 /** 00193 * @return \e f the flattening of the ellipsoid. This is the value used in 00194 * the constructor. 00195 **********************************************************************/ 00196 Math::real Flattening() const { return _f; } 00197 00198 /// \cond SKIP 00199 /** 00200 * <b>DEPRECATED</b> 00201 * @return \e r the inverse flattening of the ellipsoid. 00202 **********************************************************************/ 00203 Math::real InverseFlattening() const { return 1/_f; } 00204 /// \endcond 00205 00206 /** 00207 * @return \e k0 central scale for the projection. This is the value of \e 00208 * k0 used in the constructor and is the scale on the central meridian. 00209 **********************************************************************/ 00210 Math::real CentralScale() const { return _k0; } 00211 ///@} 00212 00213 /** 00214 * A global instantiation of TransverseMercator with the WGS84 ellipsoid 00215 * and the UTM scale factor. However, unlike UTM, no false easting or 00216 * northing is added. 00217 **********************************************************************/ 00218 static const TransverseMercator& UTM(); 00219 }; 00220 00221 } // namespace GeographicLib 00222 00223 #endif // GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP