TransverseMercator.hpp

Go to the documentation of this file.
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&uuml;ger's method which evaluates the projection and its
00030    * inverse in terms of a series.  See
00031    *  - L. Kr&uuml;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&uuml;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&quot; 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