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