1

2# -*- coding: utf-8 -*-

3

4 u'''A pure Python version of I{Karney}'s C{Exact Transverse Mercator} (ETM) projection.

5

6 Classes L{Etm}, L{ETMError} and L{ExactTransverseMercator}, transcoded from I{Karney}'s

7 C++ class U{TransverseMercatorExact<https://GeographicLib.SourceForge.io/C++/doc/

8 classGeographicLib_1_1TransverseMercatorExact.html>}, abbreviated as C{TMExact} below.

9

10 Class L{ExactTransverseMercator} provides C{Exact Transverse Mercator} projections while

11 instances of class L{Etm} represent ETM C{(easting, northing)} locations. See also

12 I{Karney}'s utility U{TransverseMercatorProj<https://GeographicLib.SourceForge.io/C++/doc/

13 TransverseMercatorProj.1.html>} and use C{"python[3] -m pygeodesy.etm ..."} to compared

14 the results.

15

16 Following is a copy of I{Karney}'s U{TransverseMercatorExact.hpp

17<https://GeographicLib.SourceForge.io/C++/doc/TransverseMercatorExact_8hpp_source.html>}

18 file C{Header}.

19

20 Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023) and licensed

21 under the MIT/X11 License. For more information, see the U{GeographicLib<https://

22 GeographicLib.SourceForge.io>} documentation.

23

24 The method entails using the U{Thompson Transverse Mercator<https://WikiPedia.org/

25 wiki/Transverse_Mercator_projection>} as an intermediate projection. The projections

26 from the intermediate coordinates to C{phi, lam} and C{x, y} are given by elliptic

27 functions. The inverse of these projections are found by Newton's method with a

28 suitable starting guess.

29

30 The relevant section of L.P. Lee's paper U{Conformal Projections Based On Jacobian

31 Elliptic Functions<https://DOI.org/10.3138/X687-1574-4325-WM62>} in part V, pp

32 67-101. The C++ implementation and notation closely follow Lee, with the following

33 exceptions::

34

35 Lee here Description

36

37 x/a xi Northing (unit Earth)

38

39 y/a eta Easting (unit Earth)

40

41 s/a sigma xi + i * eta

42

43 y x Easting

44

45 x y Northing

46

47 k e Eccentricity

48

49 k^2 mu Elliptic function parameter

50

51 k'^2 mv Elliptic function complementary parameter

52

53 m k Scale

54

55 zeta zeta Complex longitude = Mercator = chi in paper

56

57 s sigma Complex GK = zeta in paper

58

59 Minor alterations have been made in some of Lee's expressions in an attempt to

60 control round-off. For example, C{atanh(sin(phi))} is replaced by C{asinh(tan(phi))}

61 which maintains accuracy near C{phi = pi/2}. Such changes are noted in the code.

62'''

63# make sure int/int division yields float quotient, see .basics

64 from __future__ import division as _; del _ # PYCHOK semicolon

65

66 from pygeodesy.basics import map1, neg, neg_, _xinstanceof

67 from pygeodesy.constants import EPS, EPS02, PI_2, PI_4, _K0_UTM, \

68 _1_EPS, isnear0, isnear90, _0_0, \

69 _0_1, _0_5, _1_0, _2_0, _3_0, _4_0, \

70 _90_0, _180_0

71 from pygeodesy.datums import _ellipsoidal_datum, _WGS84

72 from pygeodesy.elliptic import _ALL_LAZY, Elliptic

73# from pygeodesy.errors import _incompatible # from .named

74 from pygeodesy.fmath import cbrt, hypot, hypot1, hypot2

75 from pygeodesy.fsums import Fsum, fsum1f_

76 from pygeodesy.interns import NN, _COMMASPACE_, _DASH_, _near_, _SPACE_, \

77 _spherical_, _usage

78 from pygeodesy.karney import _copyBit, _diff182, _fix90, _norm2, _norm180, \

79 _tand, _unsigned2

80# from pygeodesy.lazily import _ALL_LAZY # from .elliptic

81 from pygeodesy.named import callername, _incompatible, _NamedBase

82 from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple

83 from pygeodesy.props import deprecated_method, deprecated_property_RO, \

84 Property_RO, property_RO, _update_all, \

85 property_doc_

86 from pygeodesy.streprs import Fmt, fstr, pairs, unstr

87 from pygeodesy.units import Degrees, Scalar_

88 from pygeodesy.utily import atand, atan2d, sincos2

89 from pygeodesy.utm import _cmlon, _LLEB, _parseUTM5, _toBand, _toXtm8, \

90 _to7zBlldfn, Utm, UTMError

91

92 from math import asinh, atan2, degrees, radians, sinh, sqrt

93

94__all__ = _ALL_LAZY.etm

95__version__ = '23.08.24'

96

97_OVERFLOW = _1_EPS**2 # about 2e+31

98_TAYTOL = pow(EPS, 0.6)

99_TAYTOL2 = _TAYTOL * _2_0

100_TOL_10 = EPS * _0_1

101_TRIPS = 21 # C++ 10

102

103

104 def _overflow(x):

105 '''(INTERNAL) Like C{copysign0(OVERFLOW, B{x})}.

106 '''

107 return _copyBit(_OVERFLOW, x)

108

109

110 class ETMError(UTMError):

111 '''Exact Transverse Mercator (ETM) parse, projection or other

112 L{Etm} issue or L{ExactTransverseMercator} conversion failure.

113 '''

114 pass

115

116

117 class Etm(Utm):

118 '''Exact Transverse Mercator (ETM) coordinate, a sub-class of L{Utm},

119 a Universal Transverse Mercator (UTM) coordinate using the

120 L{ExactTransverseMercator} projection for highest accuracy.

121

122 @note: Conversion of (geodetic) lat- and longitudes to/from L{Etm}

123 coordinates is 3-4 times slower than to/from L{Utm}.

124

125 @see: Karney's U{Detailed Description<https://GeographicLib.SourceForge.io/

126 html/classGeographicLib_1_1TransverseMercatorExact.html#details>}.

127 '''

128 _Error = ETMError # see utm.UTMError

129 _exactTM = None

130

131 __init__ = Utm.__init__

132 '''New L{Etm} Exact Transverse Mercator coordinate, raising L{ETMError}s.

133

134 @see: L{Utm.__init__} for more information.

135

136 @example:

137

138 >>> import pygeodesy

139 >>> u = pygeodesy.Etm(31, 'N', 448251, 5411932)

140 '''

141

142 @property_doc_(''' the ETM projection (L{ExactTransverseMercator}).''')

143 def exactTM(self):

144 '''Get the ETM projection (L{ExactTransverseMercator}).

145 '''

146 if self._exactTM is None:

147 self.exactTM = self.datum.exactTM # ExactTransverseMercator(datum=self.datum)

148 return self._exactTM

149

150 @exactTM.setter # PYCHOK setter!

151 def exactTM(self, exactTM):

152 '''Set the ETM projection (L{ExactTransverseMercator}).

153

154 @raise ETMError: The B{C{exacTM}}'s datum incompatible

155 with this ETM coordinate's C{datum}.

156 '''

157 _xinstanceof(ExactTransverseMercator, exactTM=exactTM)

158

159 E = self.datum.ellipsoid

160 if E != exactTM.ellipsoid: # may be None

161 raise ETMError(repr(exactTM), txt=_incompatible(repr(E)))

162 self._exactTM = exactTM

163 self._scale0 = exactTM.k0

164

165 def parse(self, strETM, name=NN):

166 '''Parse a string to a similar L{Etm} instance.

167

168 @arg strETM: The ETM coordinate (C{str}),

169 see function L{parseETM5}.

170 @kwarg name: Optional instance name (C{str}),

171 overriding this name.

172

173 @return: The instance (L{Etm}).

174

175 @raise ETMError: Invalid B{C{strETM}}.

176

177 @see: Function L{pygeodesy.parseUPS5}, L{pygeodesy.parseUTM5}

178 and L{pygeodesy.parseUTMUPS5}.

179 '''

180 return parseETM5(strETM, datum=self.datum, Etm=self.classof,

181 name=name or self.name)

182

183 @deprecated_method

184 def parseETM(self, strETM): # PYCHOK no cover

185 '''DEPRECATED, use method L{Etm.parse}.

186 '''

187 return self.parse(strETM)

188

189 def toLatLon(self, LatLon=None, unfalse=True, **unused): # PYCHOK expected

190 '''Convert this ETM coordinate to an (ellipsoidal) geodetic point.

191

192 @kwarg LatLon: Optional, ellipsoidal class to return the geodetic

193 point (C{LatLon}) or C{None}.

194 @kwarg unfalse: Unfalse B{C{easting}} and B{C{northing}} if

195 C{falsed} (C{bool}).

196

197 @return: This ETM coordinate as (B{C{LatLon}}) or a

198 L{LatLonDatum5Tuple}C{(lat, lon, datum, gamma,

199 scale)} if B{C{LatLon}} is C{None}.

200

201 @raise ETMError: This ETM coordinate's C{exacTM} and this C{datum}

202 incompatible or no convergence transforming to

203 lat- and longitude.

204

205 @raise TypeError: Invalid or non-ellipsoidal B{C{LatLon}}.

206

207 @example:

208

209 >>> from pygeodesy import ellipsoidalVincenty as eV, Etm

210 >>> u = Etm(31, 'N', 448251.795, 5411932.678)

211 >>> ll = u.toLatLon(eV.LatLon) # 48°51′29.52″N, 002°17′40.20″E

212 '''

213 if not self._latlon or self._latlon._toLLEB_args != (unfalse, self.exactTM):

214 self._toLLEB(unfalse=unfalse)

215 return self._latlon5(LatLon)

216

217 def _toLLEB(self, unfalse=True, **unused): # PYCHOK signature

218 '''(INTERNAL) Compute (ellipsoidal) lat- and longitude.

219 '''

220 xTM, d = self.exactTM, self.datum

221 # double check that this and exactTM's ellipsoid match

222 if xTM._E != d.ellipsoid: # PYCHOK no cover

223 t = repr(d.ellipsoid)

224 raise ETMError(repr(xTM._E), txt=_incompatible(t))

225

226 e, n = self.eastingnorthing2(falsed=not unfalse)

227 lon0 = _cmlon(self.zone) if bool(unfalse) == self.falsed else None

228 lat, lon, g, k = xTM.reverse(e, n, lon0=lon0)

229

230 ll = _LLEB(lat, lon, datum=d, name=self.name) # utm._LLEB

231 ll._gamma = g

232 ll._scale = k

233 self._latlon5args(ll, _toBand, unfalse, xTM)

234

235 def toUtm(self): # PYCHOK signature

236 '''Copy this ETM to a UTM coordinate.

237

238 @return: The UTM coordinate (L{Utm}).

239 '''

240 return self._xcopy2(Utm)

241

242

243 class ExactTransverseMercator(_NamedBase):

244 '''Pure Python version of Karney's C++ class U{TransverseMercatorExact

245 <https://GeographicLib.SourceForge.io/C++/doc/TransverseMercatorExact_8cpp_source.html>},

246 a numerically exact transverse Mercator projection, further referred to as C{TMExact}.

247 '''

248 _datum = None # Datum

249 _E = None # Ellipsoid

250 _extendp = False # use extended domain

251# _iteration = None # ._sigmaInv2 and ._zetaInv2

252 _k0 = _K0_UTM # central scale factor

253 _lon0 = _0_0 # central meridian

254 _mu = _0_0 # ._E.e2, 1st eccentricity squared

255 _mv = _1_0 # _1_0 - ._mu

256 _raiser = False # throw Error

257 _sigmaC = None # _sigmaInv04 case

258 _zetaC = None # _zetaInv04 case

259

260 def __init__(self, datum=_WGS84, lon0=0, k0=_K0_UTM, extendp=False, name=NN, raiser=False):

261 '''New L{ExactTransverseMercator} projection.

262

263 @kwarg datum: The I{non-spherical} datum or ellipsoid (L{Datum},

264 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

265 @kwarg lon0: Central meridian, default (C{degrees180}).

266 @kwarg k0: Central scale factor (C{float}).

267 @kwarg extendp: Use the I{extended} domain (C{bool}), I{standard} otherwise.

268 @kwarg name: Optional name for the projection (C{str}).

269 @kwarg raiser: If C{True}, throw an L{ETMError} for convergence failures (C{bool}).

270

271 @raise ETMError: Near-spherical B{C{datum}} or C{ellipsoid} or invalid B{C{lon0}}

272 or B{C{k0}}.

273

274 @see: U{Constructor TransverseMercatorExact<https://GeographicLib.SourceForge.io/

275 html/classGeographicLib_1_1TransverseMercatorExact.html>} for more details,

276 especially on B{X{extendp}}.

277

278 @note: For all 255.5K U{TMcoords.dat<https://Zenodo.org/record/32470>} tests (with

279 C{0 <= lat <= 84} and C{0 <= lon}) the maximum error is C{5.2e-08 .forward}

280 (or 52 nano-meter) easting and northing and C{3.8e-13 .reverse} (or 0.38

281 pico-degrees) lat- and longitude (with Python 3.7.3+, 2.7.16+, PyPy6 3.5.3

282 and PyPy6 2.7.13, all in 64-bit on macOS 10.13.6 High Sierra C{x86_64} and

283 12.2 Monterey C{arm64} and C{"arm64_x86_64"}).

284 '''

285 if extendp:

286 self._extendp = True

287 if name:

288 self.name = name

289 if raiser:

290 self.raiser = True

291

292 self.datum = datum # invokes ._reset

293 self.k0 = k0

294 self.lon0 = lon0

295

296 @property_doc_(''' the datum (L{Datum}).''')

297 def datum(self):

298 '''Get the datum (L{Datum}) or C{None}.

299 '''

300 return self._datum

301

302 @datum.setter # PYCHOK setter!

303 def datum(self, datum):

304 '''Set the datum and ellipsoid (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

305

306 @raise ETMError: Near-spherical B{C{datum}} or C{ellipsoid}.

307 '''

308 d = _ellipsoidal_datum(datum, name=self.name) # raiser=_datum_)

309 self._reset(d)

310 self._datum = d

311

312 @Property_RO

313 def _e(self):

314 '''(INTERNAL) Get and cache C{_e}.

315 '''

316 return self._E.e

317

318 @Property_RO

319 def _1_e_90(self): # PYCHOK no cover

320 '''(INTERNAL) Get and cache C{(1 - _e) * 90}.

321 '''

322 return (_1_0 - self._e) * _90_0

323

324 @property_RO

325 def ellipsoid(self):

326 '''Get the ellipsoid (L{Ellipsoid}).

327 '''

328 return self._E

329

330 @Property_RO

331 def _e_PI_2(self):

332 '''(INTERNAL) Get and cache C{_e * PI / 2}.

333 '''

334 return self._e * PI_2

335

336 @Property_RO

337 def _e_PI_4_(self):

338 '''(INTERNAL) Get and cache C{- _e * PI / 4}.

339 '''

340 return -self._e * PI_4

341

342 @Property_RO

343 def _1_e_PI_2(self):

344 '''(INTERNAL) Get and cache C{(1 - _e) * PI / 2}.

345 '''

346 return (_1_0 - self._e) * PI_2

347

348 @Property_RO

349 def _1_2e_PI_2(self):

350 '''(INTERNAL) Get and cache C{(1 - 2 * _e) * PI / 2}.

351 '''

352 return (_1_0 - self._e * _2_0) * PI_2

353

354 @property_RO

355 def equatoradius(self):

356 '''Get this C{ellipsoid}'s equatorial radius, semi-axis (C{meter}).

357 '''

358 return self._E.a

359

360 a = equatoradius

361

362 @Property_RO

363 def _e_TAYTOL(self):

364 '''(INTERNAL) Get and cache C{e * TAYTOL}.

365 '''

366 return self._e * _TAYTOL

367

368 @Property_RO

369 def _Eu(self):

370 '''(INTERNAL) Get and cache C{Elliptic(_mu)}.

371 '''

372 return Elliptic(self._mu)

373

374 @Property_RO

375 def _Eu_cE(self):

376 '''(INTERNAL) Get and cache C{_Eu.cE}.

377 '''

378 return self._Eu.cE

379

380 def _Eu_2cE_(self, xi):

381 '''(INTERNAL) Return C{_Eu.cE * 2 - B{xi}}.

382 '''

383 return self._Eu_cE * _2_0 - xi

384

385 @Property_RO

386 def _Eu_cE_4(self):

387 '''(INTERNAL) Get and cache C{_Eu.cE / 4}.

388 '''

389 return self._Eu_cE / _4_0

390

391 @Property_RO

392 def _Eu_cK(self):

393 '''(INTERNAL) Get and cache C{_Eu.cK}.

394 '''

395 return self._Eu.cK

396

397 @Property_RO

398 def _Eu_cK_cE(self):

399 '''(INTERNAL) Get and cache C{_Eu.cK / _Eu.cE}.

400 '''

401 return self._Eu_cK / self._Eu_cE

402

403 @Property_RO

404 def _Eu_2cK_PI(self):

405 '''(INTERNAL) Get and cache C{_Eu.cK * 2 / PI}.

406 '''

407 return self._Eu_cK / PI_2

408

409 @Property_RO

410 def _Ev(self):

411 '''(INTERNAL) Get and cache C{Elliptic(_mv)}.

412 '''

413 return Elliptic(self._mv)

414

415 @Property_RO

416 def _Ev_cK(self):

417 '''(INTERNAL) Get and cache C{_Ev.cK}.

418 '''

419 return self._Ev.cK

420

421 @Property_RO

422 def _Ev_cKE(self):

423 '''(INTERNAL) Get and cache C{_Ev.cKE}.

424 '''

425 return self._Ev.cKE

426

427 @Property_RO

428 def _Ev_3cKE_4(self):

429 '''(INTERNAL) Get and cache C{_Ev.cKE * 3 / 4}.

430 '''

431 return self._Ev_cKE * 0.75

432

433 @Property_RO

434 def _Ev_5cKE_4(self):

435 '''(INTERNAL) Get and cache C{_Ev.cKE * 5 / 4}.

436 '''

437 return self._Ev_cKE * 1.25

438

439 @Property_RO

440 def extendp(self):

441 '''Get the domain (C{bool}), I{extended} or I{standard}.

442 '''

443 return self._extendp

444

445 @property_RO

446 def flattening(self):

447 '''Get this C{ellipsoid}'s flattening (C{float}).

448 '''

449 return self._E.f

450

451 f = flattening

452

453 def forward(self, lat, lon, lon0=None, name=NN): # MCCABE 13

454 '''Forward projection, from geographic to transverse Mercator.

455

456 @arg lat: Latitude of point (C{degrees}).

457 @arg lon: Longitude of point (C{degrees}).

458 @kwarg lon0: Central meridian (C{degrees180}), overriding

459 the default if not C{None}.

460 @kwarg name: Optional name (C{str}).

461

462 @return: L{Forward4Tuple}C{(easting, northing, gamma, scale)}.

463

464 @see: C{void TMExact::Forward(real lon0, real lat, real lon,

465 real &x, real &y,

466 real &gamma, real &k)}.

467

468 @raise ETMError: No convergence, thrown iff property

469 C{B{raiser}=True}.

470 '''

471 lat = _fix90(lat)

472 lon, _ = _diff182((self.lon0 if lon0 is None else lon0), lon)

473 if self.extendp:

474 backside = _lat = _lon = False

475 else: # enforce the parity

476 lat, _lat = _unsigned2(lat)

477 lon, _lon = _unsigned2(lon)

478 backside = lon > 90

479 if backside: # PYCHOK no cover

480 lon = _180_0 - lon

481 if lat == 0:

482 _lat = True

483

484 # u, v = coordinates for the Thompson TM, Lee 54

485 if lat == 90: # isnear90(lat)

486 u = self._Eu_cK

487 v = self._iteration = self._zetaC = 0

488 elif lat == 0 and lon == self._1_e_90: # PYCHOK no cover

489 u = self._iteration = self._zetaC = 0

490 v = self._Ev_cK

491 else: # tau = tan(phi), taup = sinh(psi)

492 tau, lam = _tand(lat), radians(lon)

493 u, v = self._zetaInv2(self._E.es_taupf(tau), lam)

494

495 sncndn6 = self._sncndn6(u, v)

496 y, x, _ = self._sigma3(v, *sncndn6)

497 g, k = (lon, self.k0) if isnear90(lat) else \

498 self._zetaScaled(sncndn6, ll=False)

499

500 if backside:

501 y, g = self._Eu_2cE_(y), (_180_0 - g)

502 y *= self._k0_a

503 x *= self._k0_a

504 if _lat:

505 y, g = neg_(y, g)

506 if _lon:

507 x, g = neg_(x, g)

508 return Forward4Tuple(x, y, g, k, iteration=self._iteration,

509 name=name or self.name)

510

511 def _Inv03(self, psi, dlam, _3_mv_e): # (xi, deta, _3_mv)

512 '''(INTERNAL) Partial C{_zetaInv04} or C{_sigmaInv04}, Case 2

513 '''

514 # atan2(dlam-psi, psi+dlam) + 45d gives arg(zeta - zeta0) in

515 # range [-135, 225). Subtracting 180 (multiplier is negative)

516 # makes range [-315, 45). Multiplying by 1/3 (for cube root)

517 # gives range [-105, 15). In particular the range [-90, 180]

518 # in zeta space maps to [-90, 0] in w space as required.

519 a = atan2(dlam - psi, psi + dlam) / _3_0 - PI_4

520 s, c = sincos2(a)

521 h = hypot(psi, dlam)

522 r = cbrt(h * _3_mv_e)

523 u = r * c

524 v = r * s + self._Ev_cK

525 # Error using this guess is about 0.068 * rad^(5/3)

526 return u, v, h

527

528 @property_RO

529 def iteration(self):

530 '''Get the most recent C{ExactTransverseMercator.forward}

531 or C{ExactTransverseMercator.reverse} iteration number

532 (C{int}) or C{None} if not available/applicable.

533 '''

534 return self._iteration

535

536 @property_doc_(''' the central scale factor (C{float}).''')

537 def k0(self):

538 '''Get the central scale factor (C{float}), aka I{C{scale0}}.

539 '''

540 return self._k0 # aka scale0

541

542 @k0.setter # PYCHOK setter!

543 def k0(self, k0):

544 '''Set the central scale factor (C{float}), aka I{C{scale0}}.

545

546 @raise ETMError: Invalid B{C{k0}}.

547 '''

548 k0 = Scalar_(k0=k0, Error=ETMError, low=_TOL_10, high=_1_0)

549 if self._k0 != k0:

550 ExactTransverseMercator._k0_a._update(self) # redo ._k0_a

551 self._k0 = k0

552

553 @Property_RO

554 def _k0_a(self):

555 '''(INTERNAL) Get and cache C{k0 * equatoradius}.

556 '''

557 return self.k0 * self.equatoradius

558

559 @property_doc_(''' the central meridian (C{degrees180}).''')

560 def lon0(self):

561 '''Get the central meridian (C{degrees180}).

562 '''

563 return self._lon0

564

565 @lon0.setter # PYCHOK setter!

566 def lon0(self, lon0):

567 '''Set the central meridian (C{degrees180}).

568

569 @raise ETMError: Invalid B{C{lon0}}.

570 '''

571 self._lon0 = _norm180(Degrees(lon0=lon0, Error=ETMError))

572

573 @deprecated_property_RO

574 def majoradius(self): # PYCHOK no cover

575 '''DEPRECATED, use property C{equatoradius}.'''

576 return self.equatoradius

577

578 @Property_RO

579 def _1_mu_2(self):

580 '''(INTERNAL) Get and cache C{_mu / 2 + 1}.

581 '''

582 return _1_0 + self._mu * _0_5

583

584 @Property_RO

585 def _3_mv(self):

586 '''(INTERNAL) Get and cache C{3 / _mv}.

587 '''

588 return _3_0 / self._mv

589

590 @Property_RO

591 def _3_mv_e(self):

592 '''(INTERNAL) Get and cache C{3 / (_mv * _e)}.

593 '''

594 return _3_0 / (self._mv * self._e)

595

596 def _Newton2(self, taup, lam, u, v, C, *psi): # or (xi, eta, u, v)

597 '''(INTERNAL) Invert C{_zetaInv2} or C{_sigmaInv2} using Newton's method.

598

599 @return: 2-Tuple C{(u, v)}.

600

601 @raise ETMError: No convergence.

602 '''

603 sca1, tol2 = _1_0, _TOL_10

604 if psi: # _zetaInv2

605 sca1 = sca1 / hypot1(taup) # /= chokes PyChecker

606 tol2 = tol2 / max(psi[0], _1_0)**2

607

608 _zeta3 = self._zeta3

609 _zetaDwd2 = self._zetaDwd2

610 else: # _sigmaInv2

611 _zeta3 = self._sigma3

612 _zetaDwd2 = self._sigmaDwd2

613

614 d2, r = tol2, self.raiser

615 _U_2_ = Fsum(u).fsum2_

616 _V_2_ = Fsum(v).fsum2_

617 # min iterations 2, max 6 or 7, mean 3.9 or 4.0

618 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC

619 sncndn6 = self._sncndn6(u, v)

620 du, dv = _zetaDwd2(*sncndn6)

621 T, L, _ = _zeta3(v, *sncndn6)

622 T = (taup - T) * sca1

623 L -= lam

624 u, dU = _U_2_(T * du, L * dv)

625 v, dV = _V_2_(T * dv, -L * du)

626 if d2 < tol2:

627 r = False

628 break

629 d2 = hypot2(dU, dV)

630

631 self._iteration = i

632 if r: # PYCHOK no cover

633 i = callername(up=2, underOK=True)

634 t = unstr(i, taup, lam, u, v, C=C)

635 raise ETMError(Fmt.no_convergence(d2, tol2), txt=t)

636 return u, v

637

638 @property_doc_(''' raise an L{ETMError} for convergence failures (C{bool}).''')

639 def raiser(self):

640 '''Get the error setting (C{bool}).

641 '''

642 return self._raiser

643

644 @raiser.setter # PYCHOK setter!

645 def raiser(self, raiser):

646 '''Set the error setting (C{bool}), if C{True} throw an L{ETMError}

647 for convergence failures.

648 '''

649 self._raiser = bool(raiser)

650

651 def _reset(self, datum):

652 '''(INTERNAL) Set the ellipsoid and elliptic moduli.

653

654 @arg datum: Ellipsoidal datum (C{Datum}).

655

656 @raise ETMError: Near-spherical B{C{datum}} or C{ellipsoid}.

657 '''

658 E = datum.ellipsoid

659 mu = E.e2 # .eccentricity1st2

660 mv = E.e21 # _1_0 - mu

661 if isnear0(E.e) or isnear0(mu, eps0=EPS02) \

662 or isnear0(mv, eps0=EPS02): # or sqrt(mu) != E.e

663 raise ETMError(ellipsoid=E, txt=_near_(_spherical_))

664

665 if self._datum or self._E:

666 _i = ExactTransverseMercator.iteration._uname

667 _update_all(self, _i, '_sigmaC', '_zetaC') # _under

668

669 self._E = E

670 self._mu = mu

671 self._mv = mv

672

673 def reverse(self, x, y, lon0=None, name=NN):

674 '''Reverse projection, from Transverse Mercator to geographic.

675

676 @arg x: Easting of point (C{meters}).

677 @arg y: Northing of point (C{meters}).

678 @kwarg lon0: Central meridian (C{degrees180}), overriding

679 the default if not C{None}.

680 @kwarg name: Optional name (C{str}).

681

682 @return: L{Reverse4Tuple}C{(lat, lon, gamma, scale)}.

683

684 @see: C{void TMExact::Reverse(real lon0, real x, real y,

685 real &lat, real &lon,

686 real &gamma, real &k)}

687

688 @raise ETMError: No convergence, thrown iff property

689 C{B{raiser}=True}.

690 '''

691 # undoes the steps in .forward.

692 xi = y / self._k0_a

693 eta = x / self._k0_a

694 if self.extendp:

695 backside = _lat = _lon = False

696 else: # enforce the parity

697 eta, _lon = _unsigned2(eta)

698 xi, _lat = _unsigned2(xi)

699 backside = xi > self._Eu_cE

700 if backside: # PYCHOK no cover

701 xi = self._Eu_2cE_(xi)

702

703 # u, v = coordinates for the Thompson TM, Lee 54

704 if xi or eta != self._Ev_cKE:

705 u, v = self._sigmaInv2(xi, eta)

706 else: # PYCHOK no cover

707 u = self._iteration = self._sigmaC = 0

708 v = self._Ev_cK

709

710 if v or u != self._Eu_cK:

711 g, k, lat, lon = self._zetaScaled(self._sncndn6(u, v))

712 else: # PYCHOK no cover

713 g, k, lat, lon = _0_0, self.k0, _90_0, _0_0

714

715 if backside: # PYCHOK no cover

716 lon, g = (_180_0 - lon), (_180_0 - g)

717 if _lat:

718 lat, g = neg_(lat, g)

719 if _lon:

720 lon, g = neg_(lon, g)

721 lon += self.lon0 if lon0 is None else _norm180(lon0)

722 return Reverse4Tuple(lat, _norm180(lon), g, k, # _norm180(lat)

723 iteration=self._iteration,

724 name=name or self.name)

725

726 def _scaled2(self, tau, d2, snu, cnu, dnu, snv, cnv, dnv):

727 '''(INTERNAL) C{scaled}.

728

729 @note: Argument B{C{d2}} is C{_mu * cnu**2 + _mv * cnv**2}

730 from C{._zeta3}.

731

732 @return: 2-Tuple C{(convergence, scale)}.

733

734 @see: C{void TMExact::Scale(real tau, real /*lam*/,

735 real snu, real cnu, real dnu,

736 real snv, real cnv, real dnv,

737 real &gamma, real &k)}.

738 '''

739 mu, mv = self._mu, self._mv

740 cnudnv = cnu * dnv

741 # Lee 55.12 -- negated for our sign convention. g gives

742 # the bearing (clockwise from true north) of grid north

743 g = atan2d(mv * cnv * snv * snu, cnudnv * dnu)

744 # Lee 55.13 with nu given by Lee 9.1 -- in sqrt change

745 # the numerator from (1 - snu^2 * dnv^2) to (_mv * snv^2

746 # + cnu^2 * dnv^2) to maintain accuracy near phi = 90

747 # and change the denomintor from (dnu^2 + dnv^2 - 1) to

748 # (_mu * cnu^2 + _mv * cnv^2) to maintain accuracy near

749 # phi = 0, lam = 90 * (1 - e). Similarly rewrite sqrt in

750 # 9.1 as _mv + _mu * c^2 instead of 1 - _mu * sin(phi)^2

751 if d2 > 0:

752 # originally: sec2 = 1 + tau**2 # sec(phi)^2

753 # d2 = (mu * cnu**2 + mv * cnv**2)

754 # q2 = (mv * snv**2 + cnudnv**2) / d2

755 # k = sqrt(mv + mu / sec2) * sqrt(sec2) * sqrt(q2)

756 # = sqrt(mv * sec2 + mu) * sqrt(q2)

757 # = sqrt(mv + mv * tau**2 + mu) * sqrt(q2)

758 k2 = fsum1f_(mu, mv, mv * tau**2)

759 q2 = (mv * snv**2 + cnudnv**2) / d2

760 k = (sqrt(k2) * sqrt(q2) * self.k0) if \

761 (k2 > 0 and q2 > 0) else _0_0

762 else:

763 k = _OVERFLOW

764 return g, k

765

766 def _sigma3(self, v, snu, cnu, dnu, snv, cnv, dnv):

767 '''(INTERNAL) C{sigma}.

768

769 @return: 3-Tuple C{(xi, eta, d2)}.

770

771 @see: C{void TMExact::sigma(real /*u*/, real snu, real cnu, real dnu,

772 real v, real snv, real cnv, real dnv,

773 real &xi, real &eta)}.

774

775 @raise ETMError: No convergence.

776 '''

777 mu = self._mu * cnu

778 mv = self._mv * cnv

779 # Lee 55.4 writing

780 # dnu^2 + dnv^2 - 1 = _mu * cnu^2 + _mv * cnv^2

781 d2 = cnu * mu + cnv * mv

782 mu *= snu * dnu

783 mv *= snv * dnv

784 if d2 > 0: # /= chokes PyChecker

785 mu = mu / d2

786 mv = mv / d2

787 else:

788 mu, mv = map1(_overflow, mu, mv)

789 xi = self._Eu.fE(snu, cnu, dnu) - mu

790 v -= self._Ev.fE(snv, cnv, dnv) - mv

791 return xi, v, d2

792

793 def _sigmaDwd2(self, snu, cnu, dnu, snv, cnv, dnv):

794 '''(INTERNAL) C{sigmaDwd}.

795

796 @return: 2-Tuple C{(du, dv)}.

797

798 @see: C{void TMExact::dwdsigma(real /*u*/, real snu, real cnu, real dnu,

799 real /*v*/, real snv, real cnv, real dnv,

800 real &du, real &dv)}.

801 '''

802 snuv = snu * snv

803 # Reciprocal of 55.9: dw / ds = dn(w)^2/_mv,

804 # expanding complex dn(w) using A+S 16.21.4

805 d = self._mv * (cnv**2 + self._mu * snuv**2)**2

806 r = cnv * dnu * dnv

807 i = cnu * snuv * self._mu

808 du = (r**2 - i**2) / d

809 dv = neg(_2_0 * i * r / d)

810 return du, dv

811

812 def _sigmaInv2(self, xi, eta):

813 '''(INTERNAL) Invert C{sigma} using Newton's method.

814

815 @return: 2-Tuple C{(u, v)}.

816

817 @see: C{void TMExact::sigmainv(real xi, real eta,

818 real &u, real &v)}.

819

820 @raise ETMError: No convergence.

821 '''

822 u, v, t, self._sigmaC = self._sigmaInv04(xi, eta)

823 if not t:

824 u, v = self._Newton2(xi, eta, u, v, self._sigmaC)

825 return u, v

826

827 def _sigmaInv04(self, xi, eta):

828 '''(INTERNAL) Starting point for C{sigmaInv}.

829

830 @return: 4-Tuple C{(u, v, trip, Case)}.

831

832 @see: C{bool TMExact::sigmainv0(real xi, real eta,

833 real &u, real &v)}.

834 '''

835 t = False

836 d = eta - self._Ev_cKE

837 if eta > self._Ev_5cKE_4 or (xi < d and xi < -self._Eu_cE_4):

838 # sigma as a simple pole at

839 # w = w0 = Eu.K() + i * Ev.K()

840 # and sigma is approximated by

841 # sigma = (Eu.E() + i * Ev.KE()) + 1 / (w - w0)

842 u, v = _norm2(xi - self._Eu_cE, -d)

843 u += self._Eu_cK

844 v += self._Ev_cK

845 C = 1

846

847 elif (eta > self._Ev_3cKE_4 and xi < self._Eu_cE_4) or d > 0:

848 # At w = w0 = i * Ev.K(), we have

849 # sigma = sigma0 = i * Ev.KE()

850 # sigma' = sigma'' = 0

851 # including the next term in the Taylor series gives:

852 # sigma = sigma0 - _mv / 3 * (w - w0)^3

853 # When inverting this, we map arg(w - w0) = [-pi/2, -pi/6]

854 # to arg(sigma - sigma0) = [-pi/2, pi/2] mapping arg =

855 # [-pi/2, -pi/6] to [-pi/2, pi/2]

856 u, v, h = self._Inv03(xi, d, self._3_mv)

857 t = h < _TAYTOL2

858 C = 2

859

860 else: # use w = sigma * Eu.K/Eu.E (correct in limit _e -> 0)

861 u = v = self._Eu_cK_cE

862 u *= xi

863 v *= eta

864 C = 3

865

866 return u, v, t, C

867

868 def _sncndn6(self, u, v):

869 '''(INTERNAL) Get 6-tuple C{(snu, cnu, dnu, snv, cnv, dnv)}.

870 '''

871 # snu, cnu, dnu = self._Eu.sncndn(u)

872 # snv, cnv, dnv = self._Ev.sncndn(v)

873 return self._Eu.sncndn(u) + self._Ev.sncndn(v)

874

875 def toStr(self, joined=_COMMASPACE_, **kwds): # PYCHOK signature

876 '''Return a C{str} representation.

877

878 @kwarg joined: Separator to join the attribute strings

879 (C{str} or C{None} or C{NN} for non-joined).

880 @kwarg kwds: Optional, overriding keyword arguments.

881 '''

882 d = dict(datum=self.datum.name, lon0=self.lon0,

883 k0=self.k0, extendp=self.extendp)

884 if self.name:

885 d.update(name=self.name)

886 t = pairs(d, **kwds)

887 return joined.join(t) if joined else t

888

889 def _zeta3(self, unused, snu, cnu, dnu, snv, cnv, dnv): # _sigma3 signature

890 '''(INTERNAL) C{zeta}.

891

892 @return: 3-Tuple C{(taup, lambda, d2)}.

893

894 @see: C{void TMExact::zeta(real /*u*/, real snu, real cnu, real dnu,

895 real /*v*/, real snv, real cnv, real dnv,

896 real &taup, real &lam)}

897 '''

898 e, cnu2, mv = self._e, cnu**2, self._mv

899 # Overflow value like atan(overflow) = pi/2

900 t1 = t2 = _overflow(snu)

901 # Lee 54.17 but write

902 # atanh(snu * dnv) = asinh(snu * dnv / sqrt(cnu^2 + _mv * snu^2 * snv^2))

903 # atanh(_e * snu / dnv) = asinh(_e * snu / sqrt(_mu * cnu^2 + _mv * cnv^2))

904 d1 = cnu2 + mv * (snu * snv)**2

905 if d1 > EPS02: # _EPSmin

906 t1 = snu * dnv / sqrt(d1)

907 else:

908 d1 = 0

909 d2 = self._mu * cnu2 + mv * cnv**2

910 if d2 > EPS02: # _EPSmin

911 t2 = sinh(e * asinh(e * snu / sqrt(d2)))

912 else:

913 d2 = 0

914 # psi = asinh(t1) - asinh(t2)

915 # taup = sinh(psi)

916 taup = t1 * hypot1(t2) - t2 * hypot1(t1)

917 lam = (atan2(dnu * snv, cnu * cnv) -

918 atan2(cnu * snv * e, dnu * cnv) * e) if d1 and d2 else _0_0

919 return taup, lam, d2

920

921 def _zetaDwd2(self, snu, cnu, dnu, snv, cnv, dnv):

922 '''(INTERNAL) C{zetaDwd}.

923

924 @return: 2-Tuple C{(du, dv)}.

925

926 @see: C{void TMExact::dwdzeta(real /*u*/, real snu, real cnu, real dnu,

927 real /*v*/, real snv, real cnv, real dnv,

928 real &du, real &dv)}.

929 '''

930 cnu2 = cnu**2 * self._mu

931 cnv2 = cnv**2

932 dnuv = dnu * dnv

933 dnuv2 = dnuv**2

934 snuv = snu * snv

935 snuv2 = snuv**2 * self._mu

936 # Lee 54.21 but write (see A+S 16.21.4)

937 # (1 - dnu^2 * snv^2) = (cnv^2 + _mu * snu^2 * snv^2)

938 d = self._mv * (cnv2 + snuv2)**2 # max(d, EPS02)?

939 du = cnu * dnuv * (cnv2 - snuv2) / d

940 dv = cnv * snuv * (cnu2 + dnuv2) / d

941 return du, neg(dv)

942

943 def _zetaInv2(self, taup, lam):

944 '''(INTERNAL) Invert C{zeta} using Newton's method.

945

946 @return: 2-Tuple C{(u, v)}.

947

948 @see: C{void TMExact::zetainv(real taup, real lam,

949 real &u, real &v)}.

950

951 @raise ETMError: No convergence.

952 '''

953 psi = asinh(taup)

954 u, v, t, self._zetaC = self._zetaInv04(psi, lam)

955 if not t:

956 u, v = self._Newton2(taup, lam, u, v, self._zetaC, psi)

957 return u, v

958

959 def _zetaInv04(self, psi, lam):

960 '''(INTERNAL) Starting point for C{zetaInv}.

961

962 @return: 4-Tuple C{(u, v, trip, Case)}.

963

964 @see: C{bool TMExact::zetainv0(real psi, real lam, # radians

965 real &u, real &v)}.

966 '''

967 if lam > self._1_2e_PI_2:

968 d = lam - self._1_e_PI_2

969 if psi < d and psi < self._e_PI_4_: # PYCHOK no cover

970 # N.B. this branch is normally *not* taken because psi < 0

971 # is converted psi > 0 by .forward. There's a log singularity

972 # at w = w0 = Eu.K() + i * Ev.K(), corresponding to the south

973 # pole, where we have, approximately

974 # psi = _e + i * pi/2 - _e * atanh(cos(i * (w - w0)/(1 + _mu/2)))

975 # Inverting this gives:

976 e = self._e # eccentricity

977 s, c = sincos2((PI_2 - lam) / e)

978 h, r = sinh(_1_0 - psi / e), self._1_mu_2

979 u = self._Eu_cK - r * asinh(s / hypot(c, h))

980 v = self._Ev_cK - r * atan2(c, h)

981 return u, v, False, 1

982

983 elif psi < self._e_PI_2:

984 # At w = w0 = i * Ev.K(), we have

985 # zeta = zeta0 = i * (1 - _e) * pi/2

986 # zeta' = zeta'' = 0

987 # including the next term in the Taylor series gives:

988 # zeta = zeta0 - (_mv * _e) / 3 * (w - w0)^3

989 # When inverting this, we map arg(w - w0) = [-90, 0]

990 # to arg(zeta - zeta0) = [-90, 180]

991 u, v, h = self._Inv03(psi, d, self._3_mv_e)

992 return u, v, (h < self._e_TAYTOL), 2

993

994 # Use spherical TM, Lee 12.6 -- writing C{atanh(sin(lam) /

995 # cosh(psi)) = asinh(sin(lam) / hypot(cos(lam), sinh(psi)))}.

996 # This takes care of the log singularity at C{zeta = Eu.K()},

997 # corresponding to the north pole.

998 s, c = sincos2(lam)

999 h, r = sinh(psi), self._Eu_2cK_PI

1000 # But scale to put 90, 0 on the right place

1001 u = r * atan2(h, c)

1002 v = r * asinh(s / hypot(h, c))

1003 return u, v, False, 3

1004

1005 def _zetaScaled(self, sncndn6, ll=True):

1006 '''(INTERNAL) Recompute (T, L) from (u, v) to improve accuracy of Scale.

1007

1008 @arg sncndn6: 6-Tuple C{(snu, cnu, dnu, snv, cnv, dnv)}.

1009

1010 @return: 2-Tuple C{(g, k)} if not C{B{ll}} else

1011 4-tuple C{(g, k, lat, lon)}.

1012 '''

1013 t, lam, d2 = self._zeta3(None, *sncndn6)

1014 tau = self._E.es_tauf(t)

1015 g_k = self._scaled2(tau, d2, *sncndn6)

1016 if ll:

1017 g_k += atand(tau), degrees(lam)

1018 return g_k # or (g, k, lat, lon)

1019

1020

1021 def parseETM5(strUTM, datum=_WGS84, Etm=Etm, falsed=True, name=NN):

1022 '''Parse a string representing a UTM coordinate, consisting

1023 of C{"zone[band] hemisphere easting northing"}.

1024

1025 @arg strUTM: A UTM coordinate (C{str}).

1026 @kwarg datum: Optional datum to use (L{Datum}, L{Ellipsoid},

1027 L{Ellipsoid2} or L{a_f2Tuple}).

1028 @kwarg Etm: Optional class to return the UTM coordinate

1029 (L{Etm}) or C{None}.

1030 @kwarg falsed: Both easting and northing are C{falsed} (C{bool}).

1031 @kwarg name: Optional B{C{Etm}} name (C{str}).

1032

1033 @return: The UTM coordinate (B{C{Etm}}) or if B{C{Etm}} is

1034 C{None}, a L{UtmUps5Tuple}C{(zone, hemipole, easting,

1035 northing, band)}. The C{hemipole} is the hemisphere

1036 C{'N'|'S'}.

1037

1038 @raise ETMError: Invalid B{C{strUTM}}.

1039

1040 @raise TypeError: Invalid or near-spherical B{C{datum}}.

1041

1042 @example:

1043

1044 >>> u = parseETM5('31 N 448251 5411932')

1045 >>> u.toRepr() # [Z:31, H:N, E:448251, N:5411932]

1046 >>> u = parseETM5('31 N 448251.8 5411932.7')

1047 >>> u.toStr() # 31 N 448252 5411933

1048 '''

1049 r = _parseUTM5(strUTM, datum, Etm, falsed, Error=ETMError, name=name)

1050 return r

1051

1052

1053 def toEtm8(latlon, lon=None, datum=None, Etm=Etm, falsed=True,

1054 name=NN, strict=True,

1055 zone=None, **cmoff):

1056 '''Convert a lat-/longitude point to an ETM coordinate.

1057

1058 @arg latlon: Latitude (C{degrees}) or an (ellipsoidal)

1059 geodetic C{LatLon} point.

1060 @kwarg lon: Optional longitude (C{degrees}) or C{None}.

1061 @kwarg datum: Optional datum for this ETM coordinate,

1062 overriding B{C{latlon}}'s datum (L{Datum},

1063 L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}).

1064 @kwarg Etm: Optional class to return the ETM coordinate

1065 (L{Etm}) or C{None}.

1066 @kwarg falsed: False both easting and northing (C{bool}).

1067 @kwarg name: Optional B{C{Utm}} name (C{str}).

1068 @kwarg strict: Restrict B{C{lat}} to UTM ranges (C{bool}).

1069 @kwarg zone: Optional UTM zone to enforce (C{int} or C{str}).

1070 @kwarg cmoff: DEPRECATED, use B{C{falsed}}. Offset longitude

1071 from the zone's central meridian (C{bool}).

1072

1073 @return: The ETM coordinate (B{C{Etm}}) or a

1074 L{UtmUps8Tuple}C{(zone, hemipole, easting, northing,

1075 band, datum, gamma, scale)} if B{C{Etm}} is C{None}

1076 or not B{C{falsed}}. The C{hemipole} is theC{'N'|'S'}

1077 hemisphere.

1078

1079 @raise ETMError: No convergence transforming to ETM east-

1080 and northing.

1081

1082 @raise ETMError: Invalid B{C{zone}} or near-spherical or

1083 incompatible B{C{datum}} or C{ellipsoid}.

1084

1085 @raise RangeError: If B{C{lat}} outside the valid UTM bands or

1086 if B{C{lat}} or B{C{lon}} outside the valid

1087 range and L{pygeodesy.rangerrors} set to C{True}.

1088

1089 @raise TypeError: Invalid or near-spherical B{C{datum}} or

1090 B{C{latlon}} not ellipsoidal.

1091

1092 @raise ValueError: The B{C{lon}} value is missing or B{C{latlon}}

1093 is invalid.

1094 '''

1095 z, B, lat, lon, d, f, name = _to7zBlldfn(latlon, lon, datum,

1096 falsed, name, zone,

1097 strict, ETMError, **cmoff)

1098 lon0 = _cmlon(z) if f else None

1099 x, y, g, k = d.exactTM.forward(lat, lon, lon0=lon0)

1100

1101 return _toXtm8(Etm, z, lat, x, y, B, d, g, k, f,

1102 name, latlon, d.exactTM, Error=ETMError)

1103

1104

1105 if __name__ == '__main__': # MCCABE 13

1106

1107 from pygeodesy.lazily import _ALL_MODS as _MODS, printf

1108 from sys import argv, exit as _exit

1109

1110 # mimick some of I{Karney}'s utility C{TransverseMercatorProj}

1111 _f = _r = _s = _t = False

1112 _as = argv[1:]

1113 while _as and _as[0].startswith(_DASH_):

1114 _a = _as.pop(0)

1115 if len(_a) < 2:

1116 _exit('%s: option %r invalid' % (_usage(*argv), _a))

1117 elif '-forward'.startswith(_a):

1118 _f, _r = True, False

1119 elif '-reverse'.startswith(_a):

1120 _f, _r = False, True

1121 elif '-series'.startswith(_a):

1122 _s, _t = True, False

1123 elif _a == '-t':

1124 _s, _t = False, True

1125 elif '-help'.startswith(_a):

1126 _exit(_usage(argv[0], '[-s | -t]',

1127 '[-f[orward] <lat> <lon>',

1128 '| -r[everse] <easting> <northing>',

1129 '| <lat> <lon>]',

1130 '| -h[elp]'))

1131 else:

1132 _exit('%s: option %r not supported' % (_usage(*argv), _a))

1133 if len(_as) > 1:

1134 f2 = map1(float, *_as[:2])

1135 else:

1136 _exit('%s ...: incomplete' % (_usage(*argv),))

1137

1138 if _s:

1139 tm = _MODS.ktm.KTransverseMercator()

1140 else:

1141 tm = ExactTransverseMercator(extendp=_t)

1142

1143 if _f:

1144 t = tm.forward(*f2)

1145 elif _r:

1146 t = tm.reverse(*f2)

1147 else:

1148 t = tm.forward(*f2)

1149 printf(fstr(t, sep=_SPACE_))

1150 t = tm.reverse(t.easting, t.northing)

1151 printf(fstr(t, sep=_SPACE_))

1152

1153# **) MIT License

1154#

1155# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

1156#

1157# Permission is hereby granted, free of charge, to any person obtaining a

1158# copy of this software and associated documentation files (the "Software"),

1159# to deal in the Software without restriction, including without limitation

1160# the rights to use, copy, modify, merge, publish, distribute, sublicense,

1161# and/or sell copies of the Software, and to permit persons to whom the

1162# Software is furnished to do so, subject to the following conditions:

1163#

1164# The above copyright notice and this permission notice shall be included

1165# in all copies or substantial portions of the Software.

1166#

1167# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

1168# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

1169# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

1170# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

1171# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

1172# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

1173# OTHER DEALINGS IN THE SOFTWARE.