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