GeographicLib  1.21
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-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&uuml;ger's method which evaluates the projection and its
00031    * inverse in terms of a series.  See
00032    *  - L. Kr&uuml;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&uuml;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&quot; 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