00001 /** 00002 * \file TransverseMercator.hpp 00003 * \brief Header for GeographicLib::TransverseMercator class 00004 * 00005 * Copyright (c) Charles Karney (2008, 2009, 2010) <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 6824 2010-04-19 14:25:10Z karney $" 00012 00013 #include "GeographicLib/Constants.hpp" 00014 00015 #if !defined(TM_TX_MAXPOW) 00016 /** 00017 * The order of the series approximation used in 00018 * GeographicLib::TransverseMercator. TM_TX_MAXPOW can be set to any integer 00019 * 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 * 00038 * Krüger's method has been extended from 4th to 6th order. The maximum 00039 * errors is 5 nm (ground distance) for all positions within 35 degrees of 00040 * the central meridian. The error in the convergence is 2e-15" and the 00041 * relative error in the scale is 6e-12%%. (See \ref tmerrors for the weasel 00042 * words.) The speed penalty in going to 6th order is only about 1%. 00043 * GeographicLib::TransverseMercatorExact is an alternative implementation of 00044 * the projection using exact formulas which yield accurate (to 8 nm) results 00045 * over the entire ellipsoid. 00046 * 00047 * The ellipsoid parameters and the central scale are set in the constructor. 00048 * The central meridian (which is a trivial shift of the longitude) is 00049 * specified as the \e lon0 argument of the TransverseMercator::Forward and 00050 * TransverseMercator::Reverse functions. The latitude of origin is taken to 00051 * be the equator. There is no provision in this class for specifying a 00052 * false easting or false northing or a different latitude of origin. 00053 * However these are can be simply included by the calling funtcion. For 00054 * example, the UTMUPS class applies the false easting and false northing for 00055 * the UTM projections. A more complicated example is the British National 00056 * Grid (<a href="http://www.spatialreference.org/ref/epsg/7405/"> 00057 * EPSG:7405</a>) which requires the use of a latitude of origin. This is accommodated 00058 * by (constants from 00059 * <a href="http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/guidea.html"> 00060 * A guide to coordinate systems in Great Britain</a>): 00061 \verbatim 00062 const double 00063 a = 6377563.396, b = 6356256.910, r = a/(a - b), // Airy 1830 ellipsoid 00064 k0 = 0.9996012717, lat0 = 49, lon0 = -2, // central scale and origin 00065 fe = 400000, fn = -100000; // false easting and northing 00066 // Set up basic projection 00067 const GeographicLib::TransverseMercator OSGB(a, r, k0); 00068 double x0, y0; 00069 { 00070 double gamma, k; 00071 // Transform origin point 00072 OSGB.Forward(lon0, lat0, lon0, x0, y0, gamma, k); 00073 x0 -= fe; y0 -= fn; // Combine result with false origin 00074 } 00075 double lat, lon, x, y, gamma, k; 00076 // Sample conversion from geodetic to OSGB grid 00077 std::cin >> lat >> lon; 00078 OSGB.Forward(lon0, lat, lon, x, y, gamma, k); 00079 x -= x0; y -= y0; 00080 std::cout << x << " " << y << "\n"; 00081 // Sample conversion from OSGB grid to geodetic 00082 std::cin >> x >> y; 00083 x += x0; y += y0; 00084 OSGB.Reverse(lon0, x, y, lat, lon, gamma, k); 00085 std::cout << lat << " " << lon << "\n"; 00086 \endverbatim 00087 * 00088 * See TransverseMercator.cpp for more information on the implementation. 00089 * 00090 * See \ref transversemercator for a discussion of this projection. 00091 **********************************************************************/ 00092 00093 class TransverseMercator { 00094 private: 00095 typedef Math::real real; 00096 static const int maxpow = TM_TX_MAXPOW; 00097 static const real tol, overflow; 00098 static const int numit = 5; 00099 const real _a, _r, _f, _k0, _e2, _e, _e2m, _c, _n; 00100 // _alp[0] and _bet[0] unused 00101 real _a1, _b1, _alp[maxpow + 1], _bet[maxpow + 1]; 00102 static inline real sq(real x) throw() { return x * x; } 00103 // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right 00104 static inline real tanx(real x) throw() { 00105 real t = std::tan(x); 00106 return x >= 0 ? (t >= 0 ? t : overflow) : (t < 0 ? t : -overflow); 00107 } 00108 // Return e * atanh(e * x) for f >= 0, else return 00109 // - sqrt(-e2) * atan( sqrt(-e2) * x) for f < 0 00110 inline real eatanhe(real x) const throw() { 00111 return _f >= 0 ? _e * Math::atanh(_e * x) : - _e * std::atan(_e * x); 00112 } 00113 public: 00114 00115 /** 00116 * Constructor for a ellipsoid radius \e a (meters), reciprocal flattening 00117 * \e r, and central scale factor \e k0. Setting \e r = 0 implies \e r = 00118 * inf or flattening = 0 (i.e., a sphere). Negative \e r indicates a 00119 * prolate spheroid. An exception is thrown if \e a or \e k0 is 00120 * not positive. 00121 **********************************************************************/ 00122 TransverseMercator(real a, real r, real k0); 00123 00124 /** 00125 * Convert from latitude \e lat (degrees) and longitude \e lon (degrees) to 00126 * transverse Mercator easting \e x (meters) and northing \e y (meters). 00127 * The central meridian of the transformation is \e lon0 (degrees). Also 00128 * return the meridian convergence \e gamma (degrees) and the scale \e k. 00129 * No false easting or northing is added. \e lat should be in the range 00130 * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360]. 00131 **********************************************************************/ 00132 void Forward(real lon0, real lat, real lon, 00133 real& x, real& y, real& gamma, real& k) const throw(); 00134 00135 /** 00136 * Convert from transverse Mercator easting \e x (meters) and northing \e y 00137 * (meters) to latitude \e lat (degrees) and longitude \e lon (degrees) . 00138 * The central meridian of the transformation is \e lon0 (degrees). Also 00139 * return the meridian convergence \e gamma (degrees) and the scale \e k. 00140 * No false easting or northing is added. \e lon0 should be in the range 00141 * [-180, 360]. The value of \e lon returned is in the range [-180, 180). 00142 **********************************************************************/ 00143 void Reverse(real lon0, real x, real y, 00144 real& lat, real& lon, real& gamma, real& k) const throw(); 00145 00146 /** 00147 * The major radius of the ellipsoid (meters). This is that value of \e a 00148 * used in the constructor. 00149 **********************************************************************/ 00150 Math::real MajorRadius() const throw() { return _a; } 00151 00152 /** 00153 * The inverse flattening of the ellipsoid. This is that value of \e r 00154 * used in the constructor. A value of 0 is returned for a sphere 00155 * (infinite inverse flattening). 00156 **********************************************************************/ 00157 Math::real InverseFlattening() const throw() { return _r; } 00158 00159 /** 00160 * The central scale for the projection. This is that value of \e k0 used 00161 * in the constructor and is the scale on the central meridian. 00162 **********************************************************************/ 00163 Math::real CentralScale() const throw() { return _k0; } 00164 00165 /** 00166 * A global instantiation of TransverseMercator with the WGS84 ellipsoid 00167 * and the UTM scale factor. However, unlike UTM, no false easting or 00168 * northing is added. 00169 **********************************************************************/ 00170 const static TransverseMercator UTM; 00171 }; 00172 00173 } // namespace GeographicLib 00174 00175 #endif