GeographicLib  1.35
TransverseMercatorExact.cpp
Go to the documentation of this file.
1 /**
2  * \file TransverseMercatorExact.cpp
3  * \brief Implementation for GeographicLib::TransverseMercatorExact class
4  *
5  * Copyright (c) Charles Karney (2008-2011) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  *
9  * The relevant section of Lee's paper is part V, pp 67--101,
10  * <a href="http://dx.doi.org/10.3138/X687-1574-4325-WM62">Conformal
11  * Projections Based On Jacobian Elliptic Functions</a>.
12  *
13  * The method entails using the Thompson Transverse Mercator as an
14  * intermediate projection. The projections from the intermediate
15  * coordinates to [\e phi, \e lam] and [\e x, \e y] are given by elliptic
16  * functions. The inverse of these projections are found by Newton's method
17  * with a suitable starting guess.
18  *
19  * This implementation and notation closely follows Lee, with the following
20  * exceptions:
21  * <center><table>
22  * <tr><th>Lee <th>here <th>Description
23  * <tr><td>x/a <td>xi <td>Northing (unit Earth)
24  * <tr><td>y/a <td>eta <td>Easting (unit Earth)
25  * <tr><td>s/a <td>sigma <td>xi + i * eta
26  * <tr><td>y <td>x <td>Easting
27  * <tr><td>x <td>y <td>Northing
28  * <tr><td>k <td>e <td>eccentricity
29  * <tr><td>k^2 <td>mu <td>elliptic function parameter
30  * <tr><td>k'^2 <td>mv <td>elliptic function complementary parameter
31  * <tr><td>m <td>k <td>scale
32  * <tr><td>zeta <td>zeta <td>complex longitude = Mercator = chi in paper
33  * <tr><td>s <td>sigma <td>complex GK = zeta in paper
34  * </table></center>
35  *
36  * Minor alterations have been made in some of Lee's expressions in an
37  * attempt to control round-off. For example atanh(sin(phi)) is replaced by
38  * asinh(tan(phi)) which maintains accuracy near phi = pi/2. Such changes
39  * are noted in the code.
40  **********************************************************************/
41 
43 
44 namespace GeographicLib {
45 
46  using namespace std;
47 
48  const Math::real TransverseMercatorExact::tol_ =
49  numeric_limits<real>::epsilon();
50  const Math::real TransverseMercatorExact::tol1_ = real(0.1) * sqrt(tol_);
51  const Math::real TransverseMercatorExact::tol2_ = real(0.1) * tol_;
52  const Math::real TransverseMercatorExact::taytol_ = pow(tol_, real(0.6));
53  // Overflow value s.t. atan(overflow_) = pi/2
54  const Math::real TransverseMercatorExact::overflow_ = 1 / Math::sq(tol_);
55 
57  bool extendp)
58  : _a(a)
59  , _f(f <= 1 ? f : 1/f)
60  , _k0(k0)
61  , _mu(_f * (2 - _f)) // e^2
62  , _mv(1 - _mu) // 1 - e^2
63  , _e(sqrt(_mu))
64  , _extendp(extendp)
65  , _Eu(_mu)
66  , _Ev(_mv)
67  {
68  if (!(Math::isfinite(_a) && _a > 0))
69  throw GeographicErr("Major radius is not positive");
70  if (!(_f > 0))
71  throw GeographicErr("Flattening is not positive");
72  if (!(_f < 1))
73  throw GeographicErr("Minor radius is not positive");
74  if (!(Math::isfinite(_k0) && _k0 > 0))
75  throw GeographicErr("Scale is not positive");
76  }
77 
79  TransverseMercatorExact::UTM(Constants::WGS84_a<real>(),
80  Constants::WGS84_f<real>(),
81  Constants::UTM_k0<real>());
82 
83  // tau = tan(phi), taup = sinh(psi)
84  Math::real TransverseMercatorExact::taup(real tau) const throw() {
85  real
86  tau1 = Math::hypot(real(1), tau),
87  sig = sinh( _e * Math::atanh(_e * tau / tau1) );
88  return Math::hypot(real(1), sig) * tau - sig * tau1;
89  }
90 
91  Math::real TransverseMercatorExact::taupinv(real taup) const throw() {
92  real
93  // See comment in TransverseMercator.cpp about the initial guess
94  tau = taup/_mv,
95  stol = tol_ * max(real(1), abs(taup));
96  // min iterations = 1, max iterations = 2; mean = 1.94
97  for (int i = 0; i < numit_; ++i) {
98  real
99  tau1 = Math::hypot(real(1), tau),
100  sig = sinh( _e * Math::atanh(_e * tau / tau1 ) ),
101  taupa = Math::hypot(real(1), sig) * tau - sig * tau1,
102  dtau = (taup - taupa) * (1 + _mv * Math::sq(tau)) /
103  ( _mv * tau1 * Math::hypot(real(1), taupa) );
104  tau += dtau;
105  if (!(abs(dtau) >= stol))
106  break;
107  }
108  return tau;
109  }
110 
111  void TransverseMercatorExact::zeta(real /*u*/, real snu, real cnu, real dnu,
112  real /*v*/, real snv, real cnv, real dnv,
113  real& taup, real& lam) const throw() {
114  // Lee 54.17 but write
115  // atanh(snu * dnv) = asinh(snu * dnv / sqrt(cnu^2 + _mv * snu^2 * snv^2))
116  // atanh(_e * snu / dnv) =
117  // asinh(_e * snu / sqrt(_mu * cnu^2 + _mv * cnv^2))
118  real
119  d1 = sqrt(Math::sq(cnu) + _mv * Math::sq(snu * snv)),
120  d2 = sqrt(_mu * Math::sq(cnu) + _mv * Math::sq(cnv)),
121  t1 = (d1 ? snu * dnv / d1 : (snu < 0 ? -overflow_ : overflow_)),
122  t2 = (d2 ? sinh( _e * Math::asinh(_e * snu / d2) ) :
123  (snu < 0 ? -overflow_ : overflow_));
124  // psi = asinh(t1) - asinh(t2)
125  // taup = sinh(psi)
126  taup = t1 * Math::hypot(real(1), t2) - t2 * Math::hypot(real(1), t1);
127  lam = (d1 != 0 && d2 != 0) ?
128  atan2(dnu * snv, cnu * cnv) - _e * atan2(_e * cnu * snv, dnu * cnv) :
129  0;
130  }
131 
132  void TransverseMercatorExact::dwdzeta(real /*u*/,
133  real snu, real cnu, real dnu,
134  real /*v*/,
135  real snv, real cnv, real dnv,
136  real& du, real& dv) const throw() {
137  // Lee 54.21 but write (1 - dnu^2 * snv^2) = (cnv^2 + _mu * snu^2 * snv^2)
138  // (see A+S 16.21.4)
139  real d = _mv * Math::sq(Math::sq(cnv) + _mu * Math::sq(snu * snv));
140  du = cnu * dnu * dnv * (Math::sq(cnv) - _mu * Math::sq(snu * snv)) / d;
141  dv = -snu * snv * cnv * (Math::sq(dnu * dnv) + _mu * Math::sq(cnu)) / d;
142  }
143 
144  // Starting point for zetainv
145  bool TransverseMercatorExact::zetainv0(real psi, real lam, real& u, real& v)
146  const throw() {
147  bool retval = false;
148  if (psi < -_e * Math::pi<real>()/4 &&
149  lam > (1 - 2 * _e) * Math::pi<real>()/2 &&
150  psi < lam - (1 - _e) * Math::pi<real>()/2) {
151  // N.B. this branch is normally not taken because psi < 0 is converted
152  // psi > 0 by Forward.
153  //
154  // There's a log singularity at w = w0 = Eu.K() + i * Ev.K(),
155  // corresponding to the south pole, where we have, approximately
156  //
157  // psi = _e + i * pi/2 - _e * atanh(cos(i * (w - w0)/(1 + _mu/2)))
158  //
159  // Inverting this gives:
160  real
161  psix = 1 - psi / _e,
162  lamx = (Math::pi<real>()/2 - lam) / _e;
163  u = Math::asinh(sin(lamx) / Math::hypot(cos(lamx), sinh(psix))) *
164  (1 + _mu/2);
165  v = atan2(cos(lamx), sinh(psix)) * (1 + _mu/2);
166  u = _Eu.K() - u;
167  v = _Ev.K() - v;
168  } else if (psi < _e * Math::pi<real>()/2 &&
169  lam > (1 - 2 * _e) * Math::pi<real>()/2) {
170  // At w = w0 = i * Ev.K(), we have
171  //
172  // zeta = zeta0 = i * (1 - _e) * pi/2
173  // zeta' = zeta'' = 0
174  //
175  // including the next term in the Taylor series gives:
176  //
177  // zeta = zeta0 - (_mv * _e) / 3 * (w - w0)^3
178  //
179  // When inverting this, we map arg(w - w0) = [-90, 0] to
180  // arg(zeta - zeta0) = [-90, 180]
181  real
182  dlam = lam - (1 - _e) * Math::pi<real>()/2,
183  rad = Math::hypot(psi, dlam),
184  // atan2(dlam-psi, psi+dlam) + 45d gives arg(zeta - zeta0) in range
185  // [-135, 225). Subtracting 180 (since multiplier is negative) makes
186  // range [-315, 45). Multiplying by 1/3 (for cube root) gives range
187  // [-105, 15). In particular the range [-90, 180] in zeta space maps
188  // to [-90, 0] in w space as required.
189  ang = atan2(dlam-psi, psi+dlam) - real(0.75) * Math::pi<real>();
190  // Error using this guess is about 0.21 * (rad/e)^(5/3)
191  retval = rad < _e * taytol_;
192  rad = Math::cbrt(3 / (_mv * _e) * rad);
193  ang /= 3;
194  u = rad * cos(ang);
195  v = rad * sin(ang) + _Ev.K();
196  } else {
197  // Use spherical TM, Lee 12.6 -- writing atanh(sin(lam) / cosh(psi)) =
198  // asinh(sin(lam) / hypot(cos(lam), sinh(psi))). This takes care of the
199  // log singularity at zeta = Eu.K() (corresponding to the north pole)
200  v = Math::asinh(sin(lam) / Math::hypot(cos(lam), sinh(psi)));
201  u = atan2(sinh(psi), cos(lam));
202  // But scale to put 90,0 on the right place
203  u *= _Eu.K() / (Math::pi<real>()/2);
204  v *= _Eu.K() / (Math::pi<real>()/2);
205  }
206  return retval;
207  }
208 
209  // Invert zeta using Newton's method
210  void TransverseMercatorExact::zetainv(real taup, real lam, real& u, real& v)
211  const throw() {
212  real
213  psi = Math::asinh(taup),
214  scal = 1/Math::hypot(real(1), taup);
215  if (zetainv0(psi, lam, u, v))
216  return;
217  real stol2 = tol2_ / Math::sq(max(psi, real(1)));
218  // min iterations = 2, max iterations = 6; mean = 4.0
219  for (int i = 0, trip = 0; i < numit_; ++i) {
220  real snu, cnu, dnu, snv, cnv, dnv;
221  _Eu.sncndn(u, snu, cnu, dnu);
222  _Ev.sncndn(v, snv, cnv, dnv);
223  real tau1, lam1, du1, dv1;
224  zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau1, lam1);
225  dwdzeta(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
226  tau1 -= taup;
227  lam1 -= lam;
228  tau1 *= scal;
229  real
230  delu = tau1 * du1 - lam1 * dv1,
231  delv = tau1 * dv1 + lam1 * du1;
232  u -= delu;
233  v -= delv;
234  if (trip)
235  break;
236  real delw2 = Math::sq(delu) + Math::sq(delv);
237  if (!(delw2 >= stol2))
238  ++trip;
239  }
240  }
241 
242  void TransverseMercatorExact::sigma(real /*u*/, real snu, real cnu, real dnu,
243  real v, real snv, real cnv, real dnv,
244  real& xi, real& eta) const throw() {
245  // Lee 55.4 writing
246  // dnu^2 + dnv^2 - 1 = _mu * cnu^2 + _mv * cnv^2
247  real d = _mu * Math::sq(cnu) + _mv * Math::sq(cnv);
248  xi = _Eu.E(snu, cnu, dnu) - _mu * snu * cnu * dnu / d;
249  eta = v - _Ev.E(snv, cnv, dnv) + _mv * snv * cnv * dnv / d;
250  }
251 
252  void TransverseMercatorExact::dwdsigma(real /*u*/,
253  real snu, real cnu, real dnu,
254  real /*v*/,
255  real snv, real cnv, real dnv,
256  real& du, real& dv) const throw() {
257  // Reciprocal of 55.9: dw/ds = dn(w)^2/_mv, expanding complex dn(w) using
258  // A+S 16.21.4
259  real d = _mv * Math::sq(Math::sq(cnv) + _mu * Math::sq(snu * snv));
260  real
261  dnr = dnu * cnv * dnv,
262  dni = - _mu * snu * cnu * snv;
263  du = (Math::sq(dnr) - Math::sq(dni)) / d;
264  dv = 2 * dnr * dni / d;
265  }
266 
267  // Starting point for sigmainv
268  bool TransverseMercatorExact::sigmainv0(real xi, real eta, real& u, real& v)
269  const throw() {
270  bool retval = false;
271  if (eta > real(1.25) * _Ev.KE() ||
272  (xi < -real(0.25) * _Eu.E() && xi < eta - _Ev.KE())) {
273  // sigma as a simple pole at w = w0 = Eu.K() + i * Ev.K() and sigma is
274  // approximated by
275  //
276  // sigma = (Eu.E() + i * Ev.KE()) + 1/(w - w0)
277  real
278  x = xi - _Eu.E(),
279  y = eta - _Ev.KE(),
280  r2 = Math::sq(x) + Math::sq(y);
281  u = _Eu.K() + x/r2;
282  v = _Ev.K() - y/r2;
283  } else if ((eta > real(0.75) * _Ev.KE() && xi < real(0.25) * _Eu.E())
284  || eta > _Ev.KE()) {
285  // At w = w0 = i * Ev.K(), we have
286  //
287  // sigma = sigma0 = i * Ev.KE()
288  // sigma' = sigma'' = 0
289  //
290  // including the next term in the Taylor series gives:
291  //
292  // sigma = sigma0 - _mv / 3 * (w - w0)^3
293  //
294  // When inverting this, we map arg(w - w0) = [-pi/2, -pi/6] to
295  // arg(sigma - sigma0) = [-pi/2, pi/2]
296  // mapping arg = [-pi/2, -pi/6] to [-pi/2, pi/2]
297  real
298  deta = eta - _Ev.KE(),
299  rad = Math::hypot(xi, deta),
300  // Map the range [-90, 180] in sigma space to [-90, 0] in w space. See
301  // discussion in zetainv0 on the cut for ang.
302  ang = atan2(deta-xi, xi+deta) - real(0.75) * Math::pi<real>();
303  // Error using this guess is about 0.068 * rad^(5/3)
304  retval = rad < 2 * taytol_;
305  rad = Math::cbrt(3 / _mv * rad);
306  ang /= 3;
307  u = rad * cos(ang);
308  v = rad * sin(ang) + _Ev.K();
309  } else {
310  // Else use w = sigma * Eu.K/Eu.E (which is correct in the limit _e -> 0)
311  u = xi * _Eu.K()/_Eu.E();
312  v = eta * _Eu.K()/_Eu.E();
313  }
314  return retval;
315  }
316 
317  // Invert sigma using Newton's method
318  void TransverseMercatorExact::sigmainv(real xi, real eta, real& u, real& v)
319  const throw() {
320  if (sigmainv0(xi, eta, u, v))
321  return;
322  // min iterations = 2, max iterations = 7; mean = 3.9
323  for (int i = 0, trip = 0; i < numit_; ++i) {
324  real snu, cnu, dnu, snv, cnv, dnv;
325  _Eu.sncndn(u, snu, cnu, dnu);
326  _Ev.sncndn(v, snv, cnv, dnv);
327  real xi1, eta1, du1, dv1;
328  sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi1, eta1);
329  dwdsigma(u, snu, cnu, dnu, v, snv, cnv, dnv, du1, dv1);
330  xi1 -= xi;
331  eta1 -= eta;
332  real
333  delu = xi1 * du1 - eta1 * dv1,
334  delv = xi1 * dv1 + eta1 * du1;
335  u -= delu;
336  v -= delv;
337  if (trip)
338  break;
339  real delw2 = Math::sq(delu) + Math::sq(delv);
340  if (!(delw2 >= tol2_))
341  ++trip;
342  }
343  }
344 
345  void TransverseMercatorExact::Scale(real tau, real /*lam*/,
346  real snu, real cnu, real dnu,
347  real snv, real cnv, real dnv,
348  real& gamma, real& k) const throw() {
349  real sec2 = 1 + Math::sq(tau); // sec(phi)^2
350  // Lee 55.12 -- negated for our sign convention. gamma gives the bearing
351  // (clockwise from true north) of grid north
352  gamma = atan2(_mv * snu * snv * cnv, cnu * dnu * dnv);
353  // Lee 55.13 with nu given by Lee 9.1 -- in sqrt change the numerator
354  // from
355  //
356  // (1 - snu^2 * dnv^2) to (_mv * snv^2 + cnu^2 * dnv^2)
357  //
358  // to maintain accuracy near phi = 90 and change the denomintor from
359  //
360  // (dnu^2 + dnv^2 - 1) to (_mu * cnu^2 + _mv * cnv^2)
361  //
362  // to maintain accuracy near phi = 0, lam = 90 * (1 - e). Similarly
363  // rewrite sqrt term in 9.1 as
364  //
365  // _mv + _mu * c^2 instead of 1 - _mu * sin(phi)^2
366  k = sqrt(_mv + _mu / sec2) * sqrt(sec2) *
367  sqrt( (_mv * Math::sq(snv) + Math::sq(cnu * dnv)) /
368  (_mu * Math::sq(cnu) + _mv * Math::sq(cnv)) );
369  }
370 
371  void TransverseMercatorExact::Forward(real lon0, real lat, real lon,
372  real& x, real& y, real& gamma, real& k)
373  const throw() {
375  // Explicitly enforce the parity
376  int
377  latsign = (!_extendp && lat < 0) ? -1 : 1,
378  lonsign = (!_extendp && lon < 0) ? -1 : 1;
379  lon *= lonsign;
380  lat *= latsign;
381  bool backside = !_extendp && lon > 90;
382  if (backside) {
383  if (lat == 0)
384  latsign = -1;
385  lon = 180 - lon;
386  }
387  real
388  phi = lat * Math::degree<real>(),
389  lam = lon * Math::degree<real>(),
390  tau = tanx(phi);
391 
392  // u,v = coordinates for the Thompson TM, Lee 54
393  real u, v;
394  if (lat == 90) {
395  u = _Eu.K();
396  v = 0;
397  } else if (lat == 0 && lon == 90 * (1 - _e)) {
398  u = 0;
399  v = _Ev.K();
400  } else
401  zetainv(taup(tau), lam, u, v);
402 
403  real snu, cnu, dnu, snv, cnv, dnv;
404  _Eu.sncndn(u, snu, cnu, dnu);
405  _Ev.sncndn(v, snv, cnv, dnv);
406 
407  real xi, eta;
408  sigma(u, snu, cnu, dnu, v, snv, cnv, dnv, xi, eta);
409  if (backside)
410  xi = 2 * _Eu.E() - xi;
411  y = xi * _a * _k0 * latsign;
412  x = eta * _a * _k0 * lonsign;
413 
414  if (lat == 90) {
415  gamma = lon;
416  k = 1;
417  } else {
418  // Recompute (tau, lam) from (u, v) to improve accuracy of Scale
419  zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
420  tau=taupinv(tau);
421  Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
422  gamma /= Math::degree<real>();
423  }
424  if (backside)
425  gamma = 180 - gamma;
426  gamma *= latsign * lonsign;
427  k *= _k0;
428  }
429 
430  void TransverseMercatorExact::Reverse(real lon0, real x, real y,
431  real& lat, real& lon,
432  real& gamma, real& k)
433  const throw() {
434  // This undoes the steps in Forward.
435  real
436  xi = y / (_a * _k0),
437  eta = x / (_a * _k0);
438  // Explicitly enforce the parity
439  int
440  latsign = !_extendp && y < 0 ? -1 : 1,
441  lonsign = !_extendp && x < 0 ? -1 : 1;
442  xi *= latsign;
443  eta *= lonsign;
444  bool backside = !_extendp && xi > _Eu.E();
445  if (backside)
446  xi = 2 * _Eu.E()- xi;
447 
448  // u,v = coordinates for the Thompson TM, Lee 54
449  real u, v;
450  if (xi == 0 && eta == _Ev.KE()) {
451  u = 0;
452  v = _Ev.K();
453  } else
454  sigmainv(xi, eta, u, v);
455 
456  real snu, cnu, dnu, snv, cnv, dnv;
457  _Eu.sncndn(u, snu, cnu, dnu);
458  _Ev.sncndn(v, snv, cnv, dnv);
459  real phi, lam, tau;
460  if (v != 0 || u != _Eu.K()) {
461  zeta(u, snu, cnu, dnu, v, snv, cnv, dnv, tau, lam);
462  tau = taupinv(tau);
463  phi = atan(tau);
464  lat = phi / Math::degree<real>();
465  lon = lam / Math::degree<real>();
466  Scale(tau, lam, snu, cnu, dnu, snv, cnv, dnv, gamma, k);
467  gamma /= Math::degree<real>();
468  } else {
469  tau = overflow_;
470  phi = Math::pi<real>()/2;
471  lat = 90;
472  lon = lam = gamma = 0;
473  k = 1;
474  }
475 
476  if (backside)
477  lon = 180 - lon;
478  lon *= lonsign;
479  lon = Math::AngNormalize(lon + Math::AngNormalize(lon0));
480  lat *= latsign;
481  if (backside)
482  y = 2 * _Eu.E() - y;
483  y *= _a * _k0 * latsign;
484  x *= _a * _k0 * lonsign;
485  if (backside)
486  gamma = 180 - gamma;
487  gamma *= latsign * lonsign;
488  k *= _k0;
489  }
490 
491 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:388
An exact implementation of the transverse Mercator projection.
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
static T cbrt(T x)
Definition: Math.hpp:340
static bool isfinite(T x)
Definition: Math.hpp:435
Header for GeographicLib::TransverseMercatorExact class.
static T atanh(T x)
Definition: Math.hpp:315
static T asinh(T x)
Definition: Math.hpp:288
static T hypot(T x, T y)
Definition: Math.hpp:165
TransverseMercatorExact(real a, real f, real k0, bool extendp=false)
static T sq(T x)
Definition: Math.hpp:153
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
static T AngDiff(T x, T y)
Definition: Math.hpp:418
Exception handling for GeographicLib.
Definition: Constants.hpp:320
static const TransverseMercatorExact UTM