TransverseMercatorExact.hpp

Go to the documentation of this file.
00001 /**
00002  * \file TransverseMercatorExact.hpp
00003  * \brief Header for GeographicLib::TransverseMercatorExact 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_TRANSVERSEMERCATOREXACT_HPP)
00011 #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP "$Id: TransverseMercatorExact.hpp 6950 2011-02-11 04:09:24Z karney $"
00012 
00013 #include "GeographicLib/Constants.hpp"
00014 #include "GeographicLib/EllipticFunction.hpp"
00015 
00016 namespace GeographicLib {
00017 
00018   /**
00019    * \brief An exact implementation of the Transverse Mercator Projection
00020    *
00021    * Implementation of the Transverse Mercator Projection given in
00022    *  - L. P. Lee,
00023    *    <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62"> Conformal
00024    *    Projections Based On Jacobian Elliptic Functions</a>, Part V of
00025    *    Conformal Projections Based on Elliptic Functions,
00026    *    (B. V. Gutsell, Toronto, 1976), 128pp.,
00027    *    ISBN: 0919870163
00028    *    (also appeared as:
00029    *    Monograph 16, Suppl. No. 1 to Canadian Cartographer, Vol 13).
00030    *  - C. F. F. Karney,
00031    *    <a href="http://dx.doi.org/10.1007/s00190-011-0445-3">
00032    *    Transverse Mercator with an accuracy of a few nanometers,</a>
00033    *    J. Geodesy (2011);
00034    *    preprint
00035    *    <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>.
00036    *
00037    * Lee's gives the correct results for forward and reverse
00038    * transformations subject to the branch cut rules (see the description of
00039    * the \e extendp argument to the constructor).  The maximum error is about 8
00040    * nm (ground distance) for the forward and reverse transformations.  The
00041    * error in the convergence is 2e-15&quot;, the relative error in the scale
00042    * is 7e-12%%.  See Sec. 3 of
00043    * <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for details.
00044    * The method is "exact" in the sense that the errors are close to the
00045    * round-off limit and that no changes are needed in the algorithms for them
00046    * to be used with reals of a higher precision.  Thus the errors using long
00047    * double (with a 64-bit fraction) are about 2000 times smaller than using
00048    * double (with a 53-bit fraction).
00049    *
00050    * This algorithm is about 4.5 times slower than the 6th-order Kr&uuml;ger
00051    * method, TransverseMercator, taking about 11 us for a combined forward and
00052    * reverse projection on a 2.66 GHz Intel machine (g++, version 4.3.0, -O3).
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 TransverseMercatorExact::Forward
00057    * and TransverseMercatorExact::Reverse functions.  The latitude of origin is
00058    * taken to be the equator.  See the documentation on TransverseMercator for
00059    * how to include a false easting, false northing, or a latitude of origin.
00060    *
00061    * See TransverseMercatorExact.cpp for more information on the
00062    * implementation.
00063    *
00064    * See \ref transversemercator for a discussion of this projection.
00065    **********************************************************************/
00066 
00067   class TransverseMercatorExact {
00068   private:
00069     typedef Math::real real;
00070     static const real tol, tol1, tol2, taytol, overflow;
00071     static const int numit = 10;
00072     const real _a, _r, _f, _k0, _mu, _mv, _e, _ep2;
00073     const bool _extendp;
00074     const EllipticFunction _Eu, _Ev;
00075     static inline real sq(real x) throw() { return x * x; }
00076     // tan(x) for x in [-pi/2, pi/2] ensuring that the sign is right
00077     static inline real tanx(real x) throw() {
00078       real t = std::tan(x);
00079       // Write the tests this way to ensure that tanx(NaN()) is NaN()
00080       return x >= 0 ? (!(t < 0) ? t : overflow) : (!(t >= 0) ? t : -overflow);
00081     }
00082 
00083     real taup(real tau) const throw();
00084     real taupinv(real taup) const throw();
00085 
00086     void zeta(real u, real snu, real cnu, real dnu,
00087               real v, real snv, real cnv, real dnv,
00088               real& taup, real& lam) const throw();
00089 
00090     void dwdzeta(real u, real snu, real cnu, real dnu,
00091                  real v, real snv, real cnv, real dnv,
00092                  real& du, real& dv) const throw();
00093 
00094     bool zetainv0(real psi, real lam, real& u, real& v) const throw();
00095     void zetainv(real taup, real lam, real& u, real& v) const throw();
00096 
00097     void sigma(real u, real snu, real cnu, real dnu,
00098                real v, real snv, real cnv, real dnv,
00099                real& xi, real& eta) const throw();
00100 
00101     void dwdsigma(real u, real snu, real cnu, real dnu,
00102                   real v, real snv, real cnv, real dnv,
00103                   real& du, real& dv) const throw();
00104 
00105     bool sigmainv0(real xi, real eta, real& u, real& v) const throw();
00106     void sigmainv(real xi, real eta, real& u, real& v) const throw();
00107 
00108     void Scale(real tau, real lam,
00109                real snu, real cnu, real dnu,
00110                real snv, real cnv, real dnv,
00111                real& gamma, real& k) const throw();
00112 
00113   public:
00114 
00115     /**
00116      * Constructor for a ellipsoid with
00117      *
00118      * @param[in] a equatorial radius (meters)
00119      * @param[in] r reciprocal flattening.
00120      * @param[in] k0 central scale factor.
00121      * @param[in] extendp use extended domain.
00122      *
00123      * The transverse Mercator projection has a branch point singularity at \e
00124      * lat = 0 and \e lon - \e lon0 = 90 (1 - \e e) or (for
00125      * TransverseMercatorExact::UTM) x = 18381 km, y = 0m.  The \e extendp
00126      * argument governs where the branch cut is placed.  With \e extendp =
00127      * false, the "standard" convention is followed, namely the cut is placed
00128      * along x > 18381 km, y = 0m.  Forward can be called with any \e lat and
00129      * \e lon then produces the transformation shown in Lee, Fig 46.  Reverse
00130      * analytically continues this in the +/- \e x direction.  As a
00131      * consequence, Reverse may map multiple points to the same geographic
00132      * location; for example, for TransverseMercatorExact::UTM, \e x =
00133      * 22051449.037349 m, \e y = -7131237.022729 m and \e x = 29735142.378357
00134      * m, \e y = 4235043.607933 m both map to \e lat = -2 deg, \e lon = 88 deg.
00135      *
00136      * With \e extendp = true, the branch cut is moved to the lower left
00137      * quadrant.  The various symmetries of the transverse Mercator projection
00138      * can be used to explore the projection on any sheet.  In this mode the
00139      * domains of \e lat, \e lon, \e x, and \e y are restricted to
00140      * - the union of
00141      *   - \e lat in [0, 90] and \e lon - \e lon0 in [0, 90]
00142      *   - \e lat in (-90, 0] and \e lon - \e lon0 in [90 (1 - \e e), 90]
00143      * - the union of
00144      *   - \e x/(\e k0 \e a) in [0, inf) and
00145      *     \e y/(\e k0 \e a) in [0, E(\e e^2)]
00146      *   - \e x/(\e k0 \e a) in [K(1 - \e e^2) - E(1 - \e e^2), inf) and
00147      *     \e y/(\e k0 \e a) in (-inf, 0]
00148      * .
00149      * See Sec. 5 of
00150      * <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a> for a full
00151      * discussion of the treatment of the branch cut.
00152      *
00153      * The method will work for all ellipsoids used in terrestial geodesy.  The
00154      * method cannot be applied directly to the case of a sphere (\e r = inf)
00155      * because some the constants characterizing this method diverge in that
00156      * limit, and in practise, \e r should be smaller than about
00157      * 1/numeric_limits< real >::%epsilon().  However, TransverseMercator
00158      * treats the sphere exactly.  An exception is thrown if either axis of the
00159      * ellipsoid or \e k0 is not positive or if \e r < 1.
00160      **********************************************************************/
00161     TransverseMercatorExact(real a, real r, real k0, bool extendp = false);
00162 
00163     /**
00164      * Forward projection, from geographic to transverse Mercator.
00165      *
00166      * @param[in] lon0 central meridian of the projection (degrees).
00167      * @param[in] lat latitude of point (degrees).
00168      * @param[in] lon longitude of point (degrees).
00169      * @param[out] x easting of point (meters).
00170      * @param[out] y northing of point (meters).
00171      * @param[out] gamma meridian convergence at point (degrees).
00172      * @param[out] k scale of projection at point.
00173      *
00174      * No false easting or northing is added. \e lat should be in the range
00175      * [-90, 90]; \e lon and \e lon0 should be in the range [-180, 360].
00176      **********************************************************************/
00177     void Forward(real lon0, real lat, real lon,
00178                  real& x, real& y, real& gamma, real& k) const throw();
00179 
00180     /**
00181      * Reverse projection, from transverse Mercator to geographic.
00182      *
00183      * @param[in] lon0 central meridian of the projection (degrees).
00184      * @param[in] x easting of point (meters).
00185      * @param[in] y northing of point (meters).
00186      * @param[out] lat latitude of point (degrees).
00187      * @param[out] lon longitude of point (degrees).
00188      * @param[out] gamma meridian convergence at point (degrees).
00189      * @param[out] k scale of projection at point.
00190      *
00191      * No false easting or northing is added.  \e lon0 should be in the range
00192      * [-180, 360].  The value of \e lon returned is in the range [-180, 180).
00193      **********************************************************************/
00194     void Reverse(real lon0, real x, real y,
00195                  real& lat, real& lon, real& gamma, real& k) const throw();
00196 
00197     /**
00198      * TransverseMercatorExact::Forward without returning the convergence and
00199      * scale.
00200      **********************************************************************/
00201     void Forward(real lon0, real lat, real lon,
00202                  real& x, real& y) const throw() {
00203       real gamma, k;
00204       Forward(lon0, lat, lon, x, y, gamma, k);
00205     }
00206 
00207     /**
00208      * TransverseMercatorExact::Reverse without returning the convergence and
00209      * scale.
00210      **********************************************************************/
00211     void Reverse(real lon0, real x, real y,
00212                  real& lat, real& lon) const throw() {
00213       real gamma, k;
00214       Reverse(lon0, x, y, lat, lon, gamma, k);
00215     }
00216 
00217     /** \name Inspector functions
00218      **********************************************************************/
00219     ///@{
00220     /**
00221      * @return \e a the equatorial radius of the ellipsoid (meters).  This is
00222      *   the value used in the constructor.
00223      **********************************************************************/
00224     Math::real MajorRadius() const throw() { return _a; }
00225 
00226     /**
00227      * @return \e r the inverse flattening of the ellipsoid.  This is the
00228      *   value used in the constructor.  A value of 0 is returned for a sphere
00229      *   (infinite inverse flattening).
00230      **********************************************************************/
00231     Math::real InverseFlattening() const throw() { return _r; }
00232 
00233     /**
00234      * @return \e k0 central scale for the projection.  This is the value of \e
00235      *   k0 used in the constructor and is the scale on the central meridian.
00236      **********************************************************************/
00237     Math::real CentralScale() const throw() { return _k0; }
00238     ///@}
00239 
00240     /**
00241      * A global instantiation of TransverseMercatorExact with the WGS84
00242      * ellipsoid and the UTM scale factor.  However, unlike UTM, no false
00243      * easting or northing is added.
00244      **********************************************************************/
00245     static const TransverseMercatorExact UTM;
00246   };
00247 
00248 } // namespace GeographicLib
00249 
00250 #endif