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