TransverseMercatorExact.cpp

Go to the documentation of this file.
00001 /**
00002  * \file TransverseMercatorExact.cpp
00003  * \brief Implementation for GeographicLib::TransverseMercatorExact 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  * The relevant section of Lee's paper is part V, pp 67&ndash;101,
00010  * <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62">Conformal
00011  * Projections Based On Jacobian Elliptic Functions</a>.
00012  *
00013  * The method entails using the Thompson Transverse Mercator as an
00014  * intermediate projection.  The projections from the intermediate
00015  * coordinates to [\e phi, \e lam] and [\e x, \e y] are given by elliptic
00016  * functions.  The inverse of these projections are found by Newton's method
00017  * with a suitable starting guess.
00018  *
00019  * This implementation and notation closely follows Lee, with the following
00020  * exceptions:
00021  * <center><table>
00022  * <tr><th>Lee    <th>here    <th>Description
00023  * <tr><td>x/a    <td>xi      <td>Northing (unit Earth)
00024  * <tr><td>y/a    <td>eta     <td>Easting (unit Earth)
00025  * <tr><td>s/a    <td>sigma   <td>xi + i * eta
00026  * <tr><td>y      <td>x       <td>Easting
00027  * <tr><td>x      <td>y       <td>Northing
00028  * <tr><td>k      <td>e       <td>eccentricity
00029  * <tr><td>k^2    <td>mu      <td>elliptic function parameter
00030  * <tr><td>k'^2   <td>mv      <td>elliptic function complementary parameter
00031  * <tr><td>m      <td>k       <td>scale
00032  * </table></center>
00033  *
00034  * Minor alterations have been made in some of Lee's expressions in an
00035  * attempt to control round-off.  For example atanh(sin(phi)) is replaced by
00036  * asinh(tan(phi)) which maintains accuracy near phi = pi/2.  Such changes
00037  * are noted in the code.
00038  **********************************************************************/
00039 
00040 #include "GeographicLib/TransverseMercatorExact.hpp"
00041 
00042 #define GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_CPP "$Id: TransverseMercatorExact.cpp 6813 2010-02-03 13:57:19Z karney $"
00043 
00044 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_CPP)
00045 RCSID_DECL(GEOGRAPHICLIB_TRANSVERSEMERCATOREXACT_HPP)
00046 
00047 namespace GeographicLib {
00048 
00049   using namespace std;
00050 
00051   const Math::real TransverseMercatorExact::tol =
00052     numeric_limits<real>::epsilon();
00053   const Math::real TransverseMercatorExact::tol1 = real(0.1) * sqrt(tol);
00054   const Math::real TransverseMercatorExact::tol2 = real(0.1) * tol;
00055   const Math::real TransverseMercatorExact::taytol = pow(tol, real(0.6));
00056   // Overflow value s.t. atan(overflow) = pi/2
00057   const Math::real TransverseMercatorExact::overflow = 1 / sq(tol);
00058 
00059   TransverseMercatorExact::TransverseMercatorExact(real a, real r, real k0,
00060                                                    bool extendp)
00061     : _a(a)
00062     , _r(r)
00063     , _f(1 / _r)
00064     , _k0(k0)
00065     , _mu(_f * (2 - _f))        // e^2
00066     , _mv(1 - _mu)              // 1 - e^2
00067     , _e(sqrt(_mu))
00068     , _ep2(_mu / _mv)           // e^2 / (1 - e^2)
00069     , _extendp(extendp)
00070     , _Eu(_mu)
00071     , _Ev(_mv)
00072   {
00073     if (!(_a > 0))
00074       throw GeographicErr("Major radius is not positive");
00075     if (!(_r > 0))
00076       throw GeographicErr("Inverse flattening is not positive");
00077     if (!(_k0 > 0))
00078       throw GeographicErr("Scale is not positive");
00079   }
00080 
00081   const TransverseMercatorExact
00082   TransverseMercatorExact::UTM(Constants::WGS84_a(), Constants::WGS84_r(),
00083                                Constants::UTM_k0());
00084 
00085   // tau = tan(phi), taup = sinh(psi)
00086   Math::real TransverseMercatorExact::taup(real tau) const throw() {
00087     real
00088       tau1 = Math::hypot(real(1), tau),
00089       sig = sinh( _e * Math::atanh(_e * tau / tau1) ),
00090       sig1 = Math::hypot(real(1), sig);
00091     return tau * sig1 - sig * tau1;
00092   }
00093 
00094   Math::real TransverseMercatorExact::taupinv(real taup) const throw() {
00095     real tau = taup;
00096     for (int i = 0; i < numit; ++i) {
00097       real
00098         tau1 = Math::hypot(real(1), tau),
00099         sig = sinh( _e * Math::atanh(_e * tau / tau1 ) ),
00100         sig1 =  Math::hypot(real(1), sig),
00101         dtau = - (sig1 * tau - sig * tau1 - taup) * (1 + _mv * sq(tau)) /
00102         ( (sig1 * tau1 - sig * tau) * _mv * tau1 );
00103       tau += dtau;
00104       if (abs(dtau) < tol * max(real(1), tau))
00105         break;
00106     }
00107     return tau;
00108   }
00109 
00110   void TransverseMercatorExact::zeta(real u, real snu, real cnu, real dnu,
00111                                      real v, real snv, real cnv, real dnv,
00112                                      real& taup, real& lam) const throw() {
00113     // Lee 54.17 but write
00114     // atanh(snu * dnv) = asinh(snu * dnv / sqrt(cnu^2 + _mv * snu^2 * snv^2))
00115     // atanh(_e * snu / dnv) =
00116     //         asinh(_e * snu / sqrt(_mu * cnu^2 + _mv * cnv^2))
00117     real
00118       d1 = sqrt(sq(cnu) + _mv * sq(snu * snv)),
00119       d2 = sqrt(_mu * sq(cnu) + _mv * sq(cnv)),
00120       t1 = (d1 ? snu * dnv / d1 : snu < 0 ? -overflow : overflow),
00121       t2 = (d2 ? sinh( _e * Math::asinh(_e * snu / d2) ) :
00122             snu < 0 ? -overflow : overflow);
00123     // psi = asinh(t1) - asinh(t2)
00124     // taup = sinh(psi)
00125     taup = t1 * Math::hypot(real(1), t2) - t2 * Math::hypot(real(1), t1);
00126     lam = (d1 != 0 && d2 != 0) ?
00127       atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
00128       0;
00129   }
00130 
00131   void TransverseMercatorExact::dwdzeta(real u, real snu, real cnu, real dnu,
00132                                         real v, real snv, real cnv, real dnv,
00133                                         real& du, real& dv) const throw() {
00134     // Lee 54.21 but write (1 - dnu^2 * snv^2) = (cnv^2 + _mu * snu^2 * snv^2)
00135     // (see A+S 16.21.4)
00136     real d = _mv * sq(sq(cnv) + _mu * sq(snu * snv));
00137     du =  cnu * dnu * dnv * (sq(cnv) - _mu * sq(snu * snv)) / d;
00138     dv = -snu * snv * cnv * (sq(dnu * dnv) + _mu * sq(cnu)) / d;
00139   }
00140 
00141   // Starting point for zetainv
00142   bool TransverseMercatorExact::zetainv0(real psi, real lam, real& u, real& v)
00143     const throw() {
00144     bool retval = false;
00145     if (psi < -_e * Constants::pi()/4 &&
00146         lam > (1 - 2 * _e) * Constants::pi()/2 &&
00147         psi < lam - (1 - _e) * Constants::pi()/2) {
00148       // N.B. this branch is normally not taken because psi < 0 is converted
00149       // psi > 0 by Forward.
00150       //
00151       // There's a log singularity at w = w0 = Eu.K() + i * Ev.K(),
00152       // corresponding to the south pole, where we have, approximately
00153       //
00154       //   psi = _e + i * pi/2 - _e * atanh(cos(i * (w - w0)/(1 + _mu/2)))
00155       //
00156       // Inverting this gives:
00157       real
00158         psix = 1 - psi / _e,
00159         lamx = (Constants::pi()/2 - lam) / _e;
00160       u = Math::asinh(sin(lamx) / Math::hypot(cos(lamx), sinh(psix))) *
00161         (1 + _mu/2);
00162       v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
00163       u = _Eu.K() - u;
00164       v = _Ev.K() - v;
00165     } else if (psi < _e * Constants::pi()/2 &&
00166                lam > (1 - 2 * _e) * Constants::pi()/2) {
00167       // At w = w0 = i * Ev.K(), we have
00168       //
00169       //     zeta = zeta0 = i * (1 - _e) * pi/2
00170       //     zeta' = zeta'' = 0
00171       //
00172       // including the next term in the Taylor series gives:
00173       //
00174       // zeta = zeta0 - (_mv * _e) / 3 * (w - w0)^3
00175       //
00176       // When inverting this, we map arg(w - w0) = [-90, 0] to
00177       // arg(zeta - zeta0) = [-90, 180]
00178       real
00179         dlam = lam - (1 - _e) * Constants::pi()/2,
00180         rad = Math::hypot(psi, dlam),
00181         // atan2(dlam-psi, psi+dlam) + 45d gives arg(zeta - zeta0) in range
00182         // [-135, 225).  Subtracting 180 (since multiplier is negative) makes
00183         // range [-315, 45).  Multiplying by 1/3 (for cube root) gives range
00184         // [-105, 15).  In particular the range [-90, 180] in zeta space maps
00185         // to [-90, 0] in w space as required.
00186         ang = atan2(dlam-psi, psi+dlam) - real(0.75) * Constants::pi();
00187       // Error using this guess is about 0.21 * (rad/e)^(5/3)
00188       retval = rad < _e * taytol;
00189       rad = Math::cbrt(3 / (_mv * _e) * rad);
00190       ang /= 3;
00191       u = rad * cos(ang);
00192       v = rad * sin(ang) + _Ev.K();
00193     } else {
00194       // Use spherical TM, Lee 12.6 -- writing atanh(sin(lam) / cosh(psi)) =
00195       // asinh(sin(lam) / hypot(cos(lam), sinh(psi))).  This takes care of the
00196       // log singularity at zeta = Eu.K() (corresponding to the north pole)
00197       v = Math::asinh(sin(lam) / Math::hypot(cos(lam), sinh(psi)));
00198       u = atan2(sinh(psi), cos(lam));
00199       // But scale to put 90,0 on the right place
00200       u *= _Eu.K() / (Constants::pi()/2);
00201       v *= _Eu.K() / (Constants::pi()/2);
00202     }
00203     return retval;
00204   }
00205 
00206   // Invert zeta using Newton's method
00207   void TransverseMercatorExact::zetainv(real taup, real lam, real& u, real& v)
00208     const throw()  {
00209     real
00210       psi = Math::asinh(taup),
00211       scal = 1/Math::hypot(real(1), taup);
00212     if (zetainv0(psi, lam, u, v))
00213       return;
00214     real stol2 = tol2 / sq(max(psi, real(1)));
00215     // min iterations = 2, max iterations = 6; mean = 4.0
00216     for (int i = 0, trip = 0; i < numit; ++i) {
00217       real snu, cnu, dnu, snv, cnv, dnv;
00218       _Eu.sncndn(u, snu, cnu, dnu);
00219       _Ev.sncndn(v, snv, cnv, dnv);
00220       real tau1, lam1, du1, dv1;
00221       zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
00222       dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
00223       tau1 -= taup;
00224       lam1 -= lam;
00225       tau1 *= scal;
00226       real
00227         delu = tau1 * du1 - lam1 * dv1,
00228         delv = tau1 * dv1 + lam1 * du1;
00229       u -= delu;
00230       v -= delv;
00231       if (trip)
00232         break;
00233       real delw2 = sq(delu) + sq(delv);
00234       if (delw2 < stol2)
00235         ++trip;
00236     }
00237   }
00238 
00239   void TransverseMercatorExact::sigma(real u, real snu, real cnu, real dnu,
00240                                       real v, real snv, real cnv, real dnv,
00241                                       real& xi, real& eta) const throw() {
00242     // Lee 55.4 writing
00243     // dnu^2 + dnv^2 - 1 = _mu * cnu^2 + _mv * cnv^2
00244     real d = _mu * sq(cnu) + _mv * sq(cnv);
00245     xi = _Eu.E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
00246     eta = v - _Ev.E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
00247   }
00248 
00249   void TransverseMercatorExact::dwdsigma(real u, real snu, real cnu, real dnu,
00250                                          real v, real snv, real cnv, real dnv,
00251                                          real& du, real& dv) const throw() {
00252     // Reciprocal of 55.9: dw/ds = dn(w)^2/_mv, expanding complex dn(w) using
00253     // A+S 16.21.4
00254     real d = _mv * sq(sq(cnv) + _mu * sq(snu * snv));
00255     real
00256       dnr = dnu * cnv * dnv,
00257       dni = - _mu * snu * cnu * snv;
00258     du = (sq(dnr) - sq(dni)) / d;
00259     dv = 2 * dnr * dni / d;
00260   }
00261 
00262   // Starting point for sigmainv
00263   bool TransverseMercatorExact::sigmainv0(real xi, real eta, real& u, real& v)
00264     const throw() {
00265     bool retval = false;
00266     if (eta > real(1.25) * _Ev.KE() ||
00267         (xi < -real(0.25) * _Eu.E() && xi < eta - _Ev.KE())) {
00268       // sigma as a simple pole at w = w0 = Eu.K() + i * Ev.K() and sigma is
00269       // approximated by
00270       //
00271       // sigma = (Eu.E() + i * Ev.KE()) + 1/(w - w0)
00272       real
00273         x = xi - _Eu.E(),
00274         y = eta - _Ev.KE(),
00275         r2 = sq(x) + sq(y);
00276       u = _Eu.K() + x/r2;
00277       v = _Ev.K() - y/r2;
00278     } else if ((eta > real(0.75) * _Ev.KE() && xi < real(0.25) * _Eu.E())
00279                || eta > _Ev.KE()) {
00280       // At w = w0 = i * Ev.K(), we have
00281       //
00282       //     sigma = sigma0 = i * Ev.KE()
00283       //     sigma' = sigma'' = 0
00284       //
00285       // including the next term in the Taylor series gives:
00286       //
00287       // sigma = sigma0 - _mv / 3 * (w - w0)^3
00288       //
00289       // When inverting this, we map arg(w - w0) = [-pi/2, -pi/6] to
00290       // arg(sigma - sigma0) = [-pi/2, pi/2]
00291       // mapping arg = [-pi/2, -pi/6] to [-pi/2, pi/2]
00292       real
00293         deta = eta - _Ev.KE(),
00294         rad = Math::hypot(xi, deta),
00295         // Map the range [-90, 180] in sigma space to [-90, 0] in w space.  See
00296         // discussion in zetainv0 on the cut for ang.
00297         ang = atan2(deta-xi, xi+deta) - real(0.75) * Constants::pi();
00298       // Error using this guess is about 0.068 * rad^(5/3)
00299       retval = rad < 2 * taytol;
00300       rad = Math::cbrt(3 / _mv * rad);
00301       ang /= 3;
00302       u = rad * cos(ang);
00303       v = rad * sin(ang) + _Ev.K();
00304     } else {
00305       // Else use w = sigma * Eu.K/Eu.E (which is correct in the limit _e -> 0)
00306       u = xi * _Eu.K()/_Eu.E();
00307       v = eta * _Eu.K()/_Eu.E();
00308     }
00309     return retval;
00310   }
00311 
00312   // Invert sigma using Newton's method
00313   void TransverseMercatorExact::sigmainv(real xi, real eta, real& u, real& v)
00314     const throw() {
00315     if (sigmainv0(xi, eta, u, v))
00316       return;
00317     // min iterations = 2, max iterations = 7; mean = 3.9
00318     for (int i = 0, trip = 0; i < numit; ++i) {
00319       real snu, cnu, dnu, snv, cnv, dnv;
00320       _Eu.sncndn(u, snu, cnu, dnu);
00321       _Ev.sncndn(v, snv, cnv, dnv);
00322       real xi1, eta1, du1, dv1;
00323       sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
00324       dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
00325       xi1 -= xi;
00326       eta1 -= eta;
00327       real
00328         delu = xi1 * du1 - eta1 * dv1,
00329         delv = xi1 * dv1 + eta1 * du1;
00330       u -= delu;
00331       v -= delv;
00332       if (trip)
00333         break;
00334       real delw2 = sq(delu) + sq(delv);
00335       if (delw2 < tol2)
00336         ++trip;
00337     }
00338   }
00339 
00340   void TransverseMercatorExact::Scale(real tau, real lam,
00341                                        real snu, real cnu, real dnu,
00342                                        real snv, real cnv, real dnv,
00343                                        real& gamma, real& k) const throw() {
00344     real sec2 = 1 + sq(tau);    // sec(phi)^2
00345     // Lee 55.12 -- negated for our sign convention.  gamma gives the bearing
00346     // (clockwise from true north) of grid north
00347     gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
00348     // Lee 55.13 with nu given by Lee 9.1 -- in sqrt change the numerator
00349     // from
00350     //
00351     //    (1 - snu^2 * dnv^2) to (_mv * snv^2 + cnu^2 * dnv^2)
00352     //
00353     // to maintain accuracy near phi = 90 and change the denomintor from
00354     //
00355     //    (dnu^2 + dnv^2 - 1) to (_mu * cnu^2 + _mv * cnv^2)
00356     //
00357     // to maintain accuracy near phi = 0, lam = 90 * (1 - e).  Similarly
00358     // rewrite sqrt term in 9.1 as
00359     //
00360     //    _mv + _mu * c^2 instead of 1 - _mu * sin(phi)^2
00361     k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
00362       sqrt( (_mv * sq(snv) + sq(cnu * dnv)) /
00363             (_mu * sq(cnu) + _mv * sq(cnv)) );
00364   }
00365 
00366   void TransverseMercatorExact::Forward(real lon0, real lat, real lon,
00367                                         real& x, real& y, real& gamma, real& k)
00368     const throw() {
00369     // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer)
00370     if (lon - lon0 > 180)
00371       lon -= lon0 - 360;
00372     else if (lon - lon0 <= -180)
00373       lon -= lon0 + 360;
00374     else
00375       lon -= lon0;
00376     // Now lon in (-180, 180]
00377     // Explicitly enforce the parity
00378     int
00379       latsign = !_extendp && lat < 0 ? -1 : 1,
00380       lonsign = !_extendp && lon < 0 ? -1 : 1;
00381     lon *= lonsign;
00382     lat *= latsign;
00383     bool backside = !_extendp && lon > 90;
00384     if (backside) {
00385       if (lat == 0)
00386         latsign = -1;
00387       lon = 180 - lon;
00388     }
00389     real
00390       phi = lat * Constants::degree(),
00391       lam = lon * Constants::degree(),
00392       tau = tanx(phi);
00393 
00394     // u,v = coordinates for the Thompson TM, Lee 54
00395     real u, v;
00396     if (lat == 90) {
00397       u = _Eu.K();
00398       v = 0;
00399     } else if (lat == 0 && lon == 90 * (1 - _e)) {
00400       u = 0;
00401       v = _Ev.K();
00402     } else
00403       zetainv(taup(tau), lam, u, v);
00404 
00405     real snu, cnu, dnu, snv, cnv, dnv;
00406     _Eu.sncndn(u, snu, cnu, dnu);
00407     _Ev.sncndn(v, snv, cnv, dnv);
00408 
00409     real xi, eta;
00410     sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
00411     if (backside)
00412       xi = 2 * _Eu.E() - xi;
00413     y = xi * _a * _k0 * latsign;
00414     x = eta * _a * _k0 * lonsign;
00415 
00416     // Recompute (tau, lam) from (u, v) to improve accuracy of Scale
00417     zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
00418     tau=taupinv(tau);
00419     Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
00420     gamma /= Constants::degree();
00421     if (backside)
00422       gamma = 180 - gamma;
00423     gamma *= latsign * lonsign;
00424     k *= _k0;
00425   }
00426 
00427   void TransverseMercatorExact::Reverse(real lon0, real x, real y,
00428                                         real& lat, real& lon,
00429                                         real& gamma, real& k)
00430     const throw() {
00431     // This undoes the steps in Forward.
00432     real
00433       xi = y / (_a * _k0),
00434       eta = x / (_a * _k0);
00435     // Explicitly enforce the parity
00436     int
00437       latsign = !_extendp && y < 0 ? -1 : 1,
00438       lonsign = !_extendp && x < 0 ? -1 : 1;
00439     xi *= latsign;
00440     eta *= lonsign;
00441     bool backside = !_extendp && xi > _Eu.E();
00442     if (backside)
00443       xi = 2 * _Eu.E()- xi;
00444 
00445     // u,v = coordinates for the Thompson TM, Lee 54
00446     real u, v;
00447     if (xi == 0 && eta == _Ev.KE()) {
00448       u = 0;
00449       v = _Ev.K();
00450     } else
00451       sigmainv(xi, eta, u, v);
00452 
00453     real snu, cnu, dnu, snv, cnv, dnv;
00454     _Eu.sncndn(u, snu, cnu, dnu);
00455     _Ev.sncndn(v, snv, cnv, dnv);
00456     real phi, lam, tau;
00457     if (v != 0 || u != _Eu.K()) {
00458       zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
00459       tau = taupinv(tau);
00460       phi = atan(tau);
00461       lat = phi / Constants::degree();
00462       lon = lam / Constants::degree();
00463     } else {
00464       tau = overflow;
00465       phi = Constants::pi()/2;
00466       lat = 90;
00467       lon = lam = 0;
00468     }
00469     Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
00470     gamma /= Constants::degree();
00471     if (backside)
00472       lon = 180 - lon;
00473     lon *= lonsign;
00474     // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer)
00475     if (lon + lon0 >= 180)
00476       lon += lon0 - 360;
00477     else if (lon + lon0 < -180)
00478       lon += lon0 + 360;
00479     else
00480       lon += lon0;
00481     lat *= latsign;
00482     if (backside)
00483       y = 2 * _Eu.E() - y;
00484     y *= _a * _k0 * latsign;
00485     x *= _a * _k0 * lonsign;
00486     if (backside)
00487       gamma = 180 - gamma;
00488     gamma *= latsign * lonsign;
00489     k *= _k0;
00490   }
00491 
00492 } // namespace GeographicLib

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1