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