GeographicLib
1.21
|
00001 /** 00002 * \file TransverseMercator.hpp 00003 * \brief Header for GeographicLib::TransverseMercator class 00004 * 00005 * Copyright (c) Charles Karney (2008-2011) <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 \ 00012 "$Id: 94bb078aa13d2d7392cee5498aae7df6e9914e4a $" 00013 00014 #include <GeographicLib/Constants.hpp> 00015 00016 #if !defined(TM_TX_MAXPOW) 00017 /** 00018 * The order of the series approximation used in TransverseMercator. 00019 * TM_TX_MAXPOW can be set to any integer in [4, 8]. 00020 **********************************************************************/ 00021 #define TM_TX_MAXPOW \ 00022 (GEOGRAPHICLIB_PREC == 1 ? 6 : (GEOGRAPHICLIB_PREC == 0 ? 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 00047 * 2e-15" and the relative error in the scale is 6e-12%%. See Sec. 4 of 00048 * <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for details. 00049 * The speed penalty in going to 6th order is only about 1%. 00050 * TransverseMercatorExact is an alternative implementation of the projection 00051 * using exact formulas which yield accurate (to 8 nm) results over the 00052 * entire ellipsoid. 00053 * 00054 * The ellipsoid parameters and the central scale are set in the constructor. 00055 * The central meridian (which is a trivial shift of the longitude) is 00056 * specified as the \e lon0 argument of the TransverseMercator::Forward and 00057 * TransverseMercator::Reverse functions. The latitude of origin is taken to 00058 * be the equator. There is no provision in this class for specifying a 00059 * false easting or false northing or a different latitude of origin. 00060 * However these are can be simply included by the calling function. For 00061 * example, the UTMUPS class applies the false easting and false northing for 00062 * the UTM projections. A more complicated example is the British National 00063 * Grid (<a href="http://www.spatialreference.org/ref/epsg/7405/"> 00064 * EPSG:7405</a>) which requires the use of a latitude of origin. This is 00065 * implemented by the GeographicLib::OSGB class. 00066 * 00067 * See TransverseMercator.cpp for more information on the implementation. 00068 * 00069 * See \ref transversemercator for a discussion of this projection. 00070 * 00071 * Example of use: 00072 * \include example-TransverseMercator.cpp 00073 * 00074 * <a href="TransverseMercatorProj.1.html">TransverseMercatorProj</a> is a 00075 * command-line utility providing access to the functionality of 00076 * TransverseMercator and TransverseMercatorExact. 00077 **********************************************************************/ 00078 00079 class GEOGRAPHIC_EXPORT TransverseMercator { 00080 private: 00081 typedef Math::real real; 00082 static const int maxpow_ = TM_TX_MAXPOW; 00083 static const real tol_; 00084 static const real overflow_; 00085 static const int numit_ = 5; 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 // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right 00090 static inline real tanx(real x) throw() { 00091 real t = std::tan(x); 00092 // Write the tests this way to ensure that tanx(NaN()) is NaN() 00093 return x >= 0 ? (!(t < 0) ? t : overflow_) : (!(t >= 0) ? t : -overflow_); 00094 } 00095 // Return e * atanh(e * x) for f >= 0, else return 00096 // - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 00097 inline real eatanhe(real x) const throw() { 00098 return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x); 00099 } 00100 public: 00101 00102 /** 00103 * Constructor for a ellipsoid with 00104 * 00105 * @param[in] a equatorial radius (meters). 00106 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00107 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening 00108 * to 1/\e f. 00109 * @param[in] k0 central scale factor. 00110 * 00111 * An exception is thrown if either of the axes of the ellipsoid or \e k0 00112 * is not positive. 00113 **********************************************************************/ 00114 TransverseMercator(real a, real f, real k0); 00115 00116 /** 00117 * Forward projection, from geographic to transverse Mercator. 00118 * 00119 * @param[in] lon0 central meridian of the projection (degrees). 00120 * @param[in] lat latitude of point (degrees). 00121 * @param[in] lon longitude of point (degrees). 00122 * @param[out] x easting of point (meters). 00123 * @param[out] y northing of point (meters). 00124 * @param[out] gamma meridian convergence at point (degrees). 00125 * @param[out] k scale of projection at point. 00126 * 00127 * No false easting or northing is added. \e lat should be in the range 00128 * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360]. 00129 **********************************************************************/ 00130 void Forward(real lon0, real lat, real lon, 00131 real& x, real& y, real& gamma, real& k) const throw(); 00132 00133 /** 00134 * Reverse projection, from transverse Mercator to geographic. 00135 * 00136 * @param[in] lon0 central meridian of the projection (degrees). 00137 * @param[in] x easting of point (meters). 00138 * @param[in] y northing of point (meters). 00139 * @param[out] lat latitude of point (degrees). 00140 * @param[out] lon longitude of point (degrees). 00141 * @param[out] gamma meridian convergence at point (degrees). 00142 * @param[out] k scale of projection at point. 00143 * 00144 * No false easting or northing is added. \e lon0 should be in the range 00145 * [-180, 360]. The value of \e lon returned is in the range [-180, 180). 00146 **********************************************************************/ 00147 void Reverse(real lon0, real x, real y, 00148 real& lat, real& lon, real& gamma, real& k) const throw(); 00149 00150 /** 00151 * TransverseMercator::Forward without returning the convergence and scale. 00152 **********************************************************************/ 00153 void Forward(real lon0, real lat, real lon, 00154 real& x, real& y) const throw() { 00155 real gamma, k; 00156 Forward(lon0, lat, lon, x, y, gamma, k); 00157 } 00158 00159 /** 00160 * TransverseMercator::Reverse without returning the convergence and scale. 00161 **********************************************************************/ 00162 void Reverse(real lon0, real x, real y, 00163 real& lat, real& lon) const throw() { 00164 real gamma, k; 00165 Reverse(lon0, x, y, lat, lon, gamma, k); 00166 } 00167 00168 /** \name Inspector functions 00169 **********************************************************************/ 00170 ///@{ 00171 /** 00172 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00173 * the value used in the constructor. 00174 **********************************************************************/ 00175 Math::real MajorRadius() const throw() { return _a; } 00176 00177 /** 00178 * @return \e f the flattening of the ellipsoid. This is the value used in 00179 * the constructor. 00180 **********************************************************************/ 00181 Math::real Flattening() const throw() { return _f; } 00182 00183 /// \cond SKIP 00184 /** 00185 * <b>DEPRECATED</b> 00186 * @return \e r the inverse flattening of the ellipsoid. 00187 **********************************************************************/ 00188 Math::real InverseFlattening() const throw() { return 1/_f; } 00189 /// \endcond 00190 00191 /** 00192 * @return \e k0 central scale for the projection. This is the value of \e 00193 * k0 used in the constructor and is the scale on the central meridian. 00194 **********************************************************************/ 00195 Math::real CentralScale() const throw() { return _k0; } 00196 ///@} 00197 00198 /** 00199 * A global instantiation of TransverseMercator with the WGS84 ellipsoid 00200 * and the UTM scale factor. However, unlike UTM, no false easting or 00201 * northing is added. 00202 **********************************************************************/ 00203 static const TransverseMercator UTM; 00204 }; 00205 00206 } // namespace GeographicLib 00207 00208 #endif // GEOGRAPHICLIB_TRANSVERSEMERCATOR_HPP