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 compare

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, _0_0, _0_1, _0_5, _1_0, _2_0, \

69 _3_0, _4_0, _90_0, isnear0, isnear90

70 from pygeodesy.datums import _ellipsoidal_datum, _WGS84, _EWGS84

71# from pygeodesy.ellipsoids import _EWGS84 # from .datums

72 from pygeodesy.elliptic import _ALL_LAZY, Elliptic

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

74# from pygeodesy.fsums import Fsum # from .fmath

75 from pygeodesy.fmath import cbrt, hypot, hypot1, hypot2, Fsum

76 from pygeodesy.interns import _COMMASPACE_, _near_, _SPACE_, _spherical_

77 from pygeodesy.karney import _K_2_4, _copyBit, _diff182, _fix90, \

78 _norm2, _norm180, _tand, _unsigned2

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

80 from pygeodesy.named import callername, _incompatible, _NamedBase

81 from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple

82 from pygeodesy.props import deprecated_method, deprecated_property_RO, \

83 Property_RO, property_RO, _update_all, \

84 property_doc_

85 from pygeodesy.streprs import Fmt, pairs, unstr

86 from pygeodesy.units import Degrees, Scalar_

87 from pygeodesy.utily import atan1d, atan2d, _loneg, sincos2

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

89 _to7zBlldfn, Utm, UTMError

90

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

92

93__all__ = _ALL_LAZY.etm

94__version__ = '24.10.21'

95

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

97_TAYTOL = pow(EPS, 0.6)

98_TAYTOL2 = _TAYTOL * _2_0

99_TOL_10 = EPS * _0_1

100_TRIPS = 21 # C++ 10

101

102

103 class ETMError(UTMError):

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

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

106 '''

107 pass

108

109

110 class Etm(Utm):

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

112 a Universal Transverse Mercator (UTM) coordinate using the

113 L{ExactTransverseMercator} projection for highest accuracy.

114

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

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

117

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

119 C++/doc/classGeographicLib_1_1TransverseMercatorExact.html#details>}.

120 '''

121 _Error = ETMError # see utm.UTMError

122 _exactTM = None

123

124 __init__ = Utm.__init__

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

126

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

128 '''

129

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

131 def exactTM(self):

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

133 '''

134 if self._exactTM is None:

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

136 return self._exactTM

137

138 @exactTM.setter # PYCHOK setter!

139 def exactTM(self, exactTM):

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

141

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

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

144 '''

145 _xinstanceof(ExactTransverseMercator, exactTM=exactTM)

146

147 E = self.datum.ellipsoid

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

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

150 self._exactTM = exactTM

151 self._scale0 = exactTM.k0

152

153 def parse(self, strETM, **name):

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

155

156 @arg strETM: The ETM coordinate (C{str}), see function L{parseETM5}.

157 @kwarg name: Optional C{B{name}=NN} (C{str}), overriding this name.

158

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

160

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

162

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

164 L{pygeodesy.parseUTMUPS5}.

165 '''

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

167 name=self._name__(name))

168

169 @deprecated_method

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

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

172 '''

173 return self.parse(strETM)

174

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

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

177

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

179 (C{LatLon}) or C{None}.

180 @kwarg unfalse: Unfalse B{C{easting}} and B{C{northing}} if C{falsed} (C{bool}).

181

182 @return: This ETM coordinate as (B{C{LatLon}}) or if C{B{LatLon} is None},

183 a L{LatLonDatum5Tuple}C{(lat, lon, datum, gamma, scale)}.

184

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

186 compatible or no convergence transforming to lat-/longitude.

187

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

189 '''

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

191 self._toLLEB(unfalse=unfalse)

192 return self._latlon5(LatLon)

193

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

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

196 '''

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

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

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

200 t = repr(d.ellipsoid)

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

202

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

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

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

206

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

208 self._latlon5args(ll, g, k, _toBand, unfalse, xTM)

209

210 def toUtm(self): # PYCHOK signature

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

212

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

214 '''

215 return self._xcopy2(Utm)

216

217

218 class ExactTransverseMercator(_NamedBase):

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

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

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

222 '''

223 _datum = _WGS84 # Datum

224 _E = _EWGS84 # Ellipsoid

225 _extendp = False # use extended domain

226# _iteration = None # ._sigmaInv2 and ._zetaInv2

227 _k0 = _K0_UTM # central scale factor

228 _lat0 = _0_0 # central parallel

229 _lon0 = _0_0 # central meridian

230 _mu = _EWGS84.e2 # 1st eccentricity squared

231 _mv = _EWGS84.e21 # 1 - ._mu

232 _raiser = False # throw Error

233 _sigmaC = None # most recent _sigmaInv04 case C{int}

234 _zetaC = None # most recent _zetaInv04 case C{int}

235

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

237 '''New L{ExactTransverseMercator} projection.

238

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

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

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

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

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

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

245 @kwarg name: Optional C{B{name}=NN} for the projection (C{str}).

246

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

248 or B{C{k0}}.

249

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

251 C++/doc/classGeographicLib_1_1TransverseMercatorExact.html>} for more details,

252 especially on B{X{extendp}}.

253

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

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

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

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

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

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

260 '''

261 if extendp:

262 self._extendp = True

263 if name:

264 self.name = name

265 if raiser:

266 self.raiser = True

267

268 TM = ExactTransverseMercator

269 if datum not in (TM._datum, TM._E, None):

270 self.datum = datum # invokes ._resets

271 if lon0 or lon0 != TM._lon0:

272 self.lon0 = lon0

273 if k0 is not TM._k0:

274 self.k0 = k0

275

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

277 def datum(self):

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

279 '''

280 return self._datum

281

282 @datum.setter # PYCHOK setter!

283 def datum(self, datum):

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

285

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

287 '''

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

289 self._resets(d)

290 self._datum = d

291

292 @Property_RO

293 def _e(self):

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

295 '''

296 return self._E.e

297

298 @Property_RO

299 def _1_e_90(self): # PYCHOK no cover

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

301 '''

302 return (_1_0 - self._e) * _90_0

303

304 @property_RO

305 def ellipsoid(self):

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

307 '''

308 return self._E

309

310 @Property_RO

311 def _e_PI_2(self):

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

313 '''

314 return self._e * PI_2

315

316 @Property_RO

317 def _e_PI_4_(self):

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

319 '''

320 return -self._e * PI_4

321

322 @Property_RO

323 def _1_e_PI_2(self):

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

325 '''

326 return (_1_0 - self._e) * PI_2

327

328 @Property_RO

329 def _1_2e_PI_2(self):

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

331 '''

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

333

334 @property_RO

335 def equatoradius(self):

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

337 '''

338 return self._E.a

339

340 a = equatoradius

341

342 @Property_RO

343 def _e_TAYTOL(self):

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

345 '''

346 return self._e * _TAYTOL

347

348 @Property_RO

349 def _Eu(self):

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

351 '''

352 return Elliptic(self._mu)

353

354 @Property_RO

355 def _Eu_cE(self):

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

357 '''

358 return self._Eu.cE

359

360 def _Eu_2cE_(self, xi):

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

362 '''

363 return self._Eu_cE * _2_0 - xi

364

365 @Property_RO

366 def _Eu_cE_4(self):

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

368 '''

369 return self._Eu_cE / _4_0

370

371 @Property_RO

372 def _Eu_cK(self):

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

374 '''

375 return self._Eu.cK

376

377 @Property_RO

378 def _Eu_cK_cE(self):

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

380 '''

381 return self._Eu_cK / self._Eu_cE

382

383 @Property_RO

384 def _Eu_2cK_PI(self):

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

386 '''

387 return self._Eu_cK / PI_2

388

389 @Property_RO

390 def _Ev(self):

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

392 '''

393 return Elliptic(self._mv)

394

395 @Property_RO

396 def _Ev_cK(self):

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

398 '''

399 return self._Ev.cK

400

401 @Property_RO

402 def _Ev_cKE(self):

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

404 '''

405 return self._Ev.cKE

406

407 @Property_RO

408 def _Ev_3cKE_4(self):

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

410 '''

411 return self._Ev_cKE * 0.75 # _0_75

412

413 @Property_RO

414 def _Ev_5cKE_4(self):

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

416 '''

417 return self._Ev_cKE * 1.25 # _1_25

418

419 @Property_RO

420 def extendp(self):

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

422 '''

423 return self._extendp

424

425 @property_RO

426 def flattening(self):

427 '''Get the C{ellipsoid}'s flattening (C{scalar}).

428 '''

429 return self._E.f

430

431 f = flattening

432

433 def forward(self, lat, lon, lon0=None, jam=_K_2_4, **name): # MCCABE 13

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

435

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

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

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

439 the default if not C{None}.

440 @kwarg jam: If C{True}, use the C{Jacobi amplitude}

441 otherwise C{Bulirsch}' function (C{bool}).

442 @kwarg name: Optional C{B{name}=NN} (C{str}).

443

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

445

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

447 real &x, real &y,

448 real &gamma, real &k)}.

449

450 @raise ETMError: No convergence, thrown iff property

451 C{B{raiser}=True}.

452 '''

453 lat = _fix90(lat - self._lat0)

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

455 if self.extendp:

456 backside = _lat = _lon = False

457 else: # enforce the parity

458 lat, _lat = _unsigned2(lat)

459 lon, _lon = _unsigned2(lon)

460 backside = lon > 90

461 if backside: # PYCHOK no cover

462 lon = _loneg(lon)

463 if lat == 0:

464 _lat = True

465

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

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

468 u = self._Eu_cK

469 v = self._iteration = self._zetaC = 0

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

471 u = self._iteration = self._zetaC = 0

472 v = self._Ev_cK

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

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

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

476

477 sncndn6 = self._sncndn6(u, v, jam=jam)

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

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

480 self._zetaScaled(sncndn6, ll=False)

481

482 if backside:

483 y, g = self._Eu_2cE_(y), _loneg(g)

484 y *= self._k0_a

485 x *= self._k0_a

486 if _lat:

487 y, g = neg_(y, g)

488 if _lon:

489 x, g = neg_(x, g)

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

491 name=self._name__(name))

492

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

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

495 '''

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

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

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

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

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

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

502 s, c = sincos2(a)

503 h = hypot(psi, dlam)

504 r = cbrt(h * _3_mv_e)

505 u = r * c

506 v = r * s + self._Ev_cK

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

508 return u, v, h

509

510 @property_RO

511 def iteration(self):

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

513 or C{ExactTransverseMercator.reverse} iteration number

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

515 '''

516 return self._iteration

517

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

519 def k0(self):

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

521 '''

522 return self._k0 # aka scale0

523

524 @k0.setter # PYCHOK setter!

525 def k0(self, k0):

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

527

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

529 '''

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

531 if self._k0 != k0:

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

533 self._k0 = k0

534

535 @Property_RO

536 def _k0_a(self):

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

538 '''

539 return self.k0 * self.equatoradius

540

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

542 def lon0(self):

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

544 '''

545 return self._lon0

546

547 @lon0.setter # PYCHOK setter!

548 def lon0(self, lon0):

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

550

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

552 '''

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

554

555 @deprecated_property_RO

556 def majoradius(self): # PYCHOK no cover

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

558 return self.equatoradius

559

560 @Property_RO

561 def _1_mu_2(self):

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

563 '''

564 return self._mu * _0_5 + _1_0

565

566 @Property_RO

567 def _3_mv(self):

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

569 '''

570 return _3_0 / self._mv

571

572 @Property_RO

573 def _3_mv_e(self):

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

575 '''

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

577

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

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

580

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

582

583 @raise ETMError: No convergence.

584 '''

585 sca1, tol2 = _1_0, _TOL_10

586 if psi: # _zetaInv2

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

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

589

590 _zeta3 = self._zeta3

591 _zetaDwd2 = self._zetaDwd2

592 else: # _sigmaInv2

593 _zeta3 = self._sigma3

594 _zetaDwd2 = self._sigmaDwd2

595

596 d2, r = tol2, self.raiser

597 _U_2 = Fsum(u).fsum2f_

598 _V_2 = Fsum(v).fsum2f_

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

600 _hy2 = hypot2

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

602 sncndn6 = self._sncndn6(u, v)

603 du, dv = _zetaDwd2(*sncndn6)

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

605 T = (taup - T) * sca1

606 L -= lam

607 u, dU = _U_2(T * du, L * dv)

608 v, dV = _V_2(T * dv, -L * du)

609 if d2 < tol2:

610 r = False

611 break

612 d2 = _hy2(dU, dV)

613

614 self._iteration = i

615 if r: # PYCHOK no cover

616 n = callername(up=2, underOK=True)

617 t = unstr(n, taup, lam, u, v, C=C)

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

619 return u, v

620

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

622 def raiser(self):

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

624 '''

625 return self._raiser

626

627 @raiser.setter # PYCHOK setter!

628 def raiser(self, raiser):

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

630 for convergence failures.

631 '''

632 self._raiser = bool(raiser)

633

634 def reset(self, lat0, lon0):

635 '''Set the central parallel and meridian.

636

637 @arg lat0: Latitude of the central parallel (C{degrees90}).

638 @arg lon0: Longitude of the central parallel (C{degrees180}).

639

640 @return: 2-Tuple C{(lat0, lon0)} of the previous central

641 parallel and meridian.

642

643 @raise ETMError: Invalid B{C{lat0}} or B{C{lon0}}.

644 '''

645 t = self._lat0, self.lon0

646 self._lat0 = _fix90(Degrees(lat0=lat0, Error=ETMError))

647 self. lon0 = lon0

648 return t

649

650 def _resets(self, datum):

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

652

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

654

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

656 '''

657 E = datum.ellipsoid

658 mu = E.e2 # .eccentricity1st2

659 mv = E.e21 # _1_0 - mu

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

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

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

663

664 if self._datum or self._E:

665 _i = ExactTransverseMercator.iteration._uname

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

667

668 self._E = E

669 self._mu = mu

670 self._mv = mv

671

672 def reverse(self, x, y, lon0=None, jam=_K_2_4, **name):

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

674

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

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

677 @kwarg lon0: Optional central meridian (C{degrees180}),

678 overriding the default (C{iff not None}).

679 @kwarg jam: If C{True}, use the C{Jacobi amplitude}

680 otherwise C{Bulirsch}' function (C{bool}).

681 @kwarg name: Optional C{B{name}=NN} (C{str}).

682

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

684

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

686 real &lat, real &lon,

687 real &gamma, real &k)}

688

689 @raise ETMError: No convergence, thrown iff property

690 C{B{raiser}=True}.

691 '''

692 # undoes the steps in .forward.

693 xi = y / self._k0_a

694 eta = x / self._k0_a

695 if self.extendp:

696 backside = _lat = _lon = False

697 else: # enforce the parity

698 eta, _lon = _unsigned2(eta)

699 xi, _lat = _unsigned2(xi)

700 backside = xi > self._Eu_cE

701 if backside: # PYCHOK no cover

702 xi = self._Eu_2cE_(xi)

703

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

705 if xi or eta != self._Ev_cKE:

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

707 else: # PYCHOK no cover

708 u = self._iteration = self._sigmaC = 0

709 v = self._Ev_cK

710

711 if v or u != self._Eu_cK:

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

713 else: # PYCHOK no cover

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

715

716 if backside: # PYCHOK no cover

717 lon, g = _loneg(lon), _loneg(g)

718 if _lat:

719 lat, g = neg_(lat, g)

720 if _lon:

721 lon, g = neg_(lon, g)

722 lat += self._lat0

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

724 return Reverse4Tuple(lat, _norm180(lon), g, k, # _fix90(lat)

725 iteration=self._iteration,

726 name=self._name__(name))

727

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

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

730

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

732 from C{._zeta3}.

733

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

735

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

737 real snu, real cnu, real dnu,

738 real snv, real cnv, real dnv,

739 real &gamma, real &k)}.

740 '''

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

742 cnudnv = cnu * dnv

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

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

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

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

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

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

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

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

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

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

753 if d2 > 0:

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

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

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

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

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

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

760 k, q2 = _0_0, (mv * snv**2 + cnudnv**2)

761 if q2 > 0:

762 k2 = (tau**2 + _1_0) * mv + mu

763 if k2 > 0:

764 k = sqrt(k2) * sqrt(q2 / d2) * self.k0

765 else:

766 k = _OVERFLOW

767 return g, k

768

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

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

771

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

773

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

775 real v, real snv, real cnv, real dnv,

776 real &xi, real &eta)}.

777

778 @raise ETMError: No convergence.

779 '''

780 mu = self._mu * cnu

781 mv = self._mv * cnv

782 # Lee 55.4 writing

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

784 d2 = cnu * mu + cnv * mv

785 mu *= snu * dnu

786 mv *= snv * dnv

787 if d2 > 0: # /= chokes PyChecker

788 mu = mu / d2

789 mv = mv / d2

790 else:

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

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

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

794 return xi, v, d2

795

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

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

798

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

800

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

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

803 real &du, real &dv)}.

804 '''

805 mu = self._mu

806 snuv = snu * snv

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

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

809 d = (cnv**2 + snuv**2 * mu)**2 * self._mv

810 r = cnv * dnu * dnv

811 i = cnu * snuv * mu

812 du = (r**2 - i**2) / d # (r + i) * (r - i) / d

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

814 return du, dv

815

816 def _sigmaInv2(self, xi, eta):

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

818

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

820

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

822 real &u, real &v)}.

823

824 @raise ETMError: No convergence.

825 '''

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

827 if not t:

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

829 return u, v

830

831 def _sigmaInv04(self, xi, eta):

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

833

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

835

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

837 real &u, real &v)}.

838 '''

839 t = False

840 d = eta - self._Ev_cKE

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

842 # sigma as a simple pole at

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

844 # and sigma is approximated by

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

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

847 u += self._Eu_cK

848 v += self._Ev_cK

849 C = 1

850

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

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

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

854 # sigma' = sigma'' = 0

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

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

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

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

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

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

861 t = h < _TAYTOL2

862 C = 2

863

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

865 u = v = self._Eu_cK_cE

866 u *= xi

867 v *= eta

868 C = 3

869

870 return u, v, t, C

871

872 def _sncndn6(self, u, v, **jam):

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

874 '''

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

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

877 return self._Eu.sncndn(u, **jam) + self._Ev.sncndn(v, **jam)

878

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

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

881

882 @kwarg joined: Separator to join the attribute strings

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

884 @kwarg kwds: Optional, overriding keyword arguments.

885 '''

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

887 k0=self.k0, extendp=self.extendp)

888 if self.name:

889 d.update(name=self.name)

890 t = pairs(d, **kwds)

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

892

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

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

895

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

897

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

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

900 real &taup, real &lam)}

901 '''

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

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

904 t1 = t2 = _overflow(snu)

905 # Lee 54.17 but write

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

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

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

909 if d1 > EPS02: # _EPSmin

910 t1 = snu * dnv / sqrt(d1)

911 else:

912 d1 = 0

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

914 if d2 > EPS02: # _EPSmin

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

916 else:

917 d2 = 0

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

919 # taup = sinh(psi)

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

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

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

923 return taup, lam, d2

924

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

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

927

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

929

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

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

932 real &du, real &dv)}.

933 '''

934 cnu2 = cnu**2 * self._mu

935 cnv2 = cnv**2

936 dnuv = dnu * dnv

937 dnuv2 = dnuv**2

938 snuv = snu * snv

939 snuv2 = snuv**2 * self._mu

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

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

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

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

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

945 return du, neg(dv)

946

947 def _zetaInv2(self, taup, lam):

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

949

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

951

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

953 real &u, real &v)}.

954

955 @raise ETMError: No convergence.

956 '''

957 psi = asinh(taup)

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

959 if not t:

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

961 return u, v

962

963 def _zetaInv04(self, psi, lam):

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

965

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

967

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

969 real &u, real &v)}.

970 '''

971 if lam > self._1_2e_PI_2:

972 d = lam - self._1_e_PI_2

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

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

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

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

977 # pole, where we have, approximately

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

979 # Inverting this gives:

980 e = self._e # eccentricity

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

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

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

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

985 return u, v, False, 1

986

987 elif psi < self._e_PI_2:

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

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

990 # zeta' = zeta'' = 0

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

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

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

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

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

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

997

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

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

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

1001 # corresponding to the north pole.

1002 s, c = sincos2(lam)

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

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

1005 u = r * atan2(h, c)

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

1007 return u, v, False, 3

1008

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

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

1011

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

1013

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

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

1016 '''

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

1018 tau = self._E.es_tauf(t)

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

1020 if ll:

1021 g_k += atan1d(tau), degrees(lam)

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

1023

1024

1025 def _overflow(x):

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

1027 '''

1028 return _copyBit(_OVERFLOW, x)

1029

1030

1031 def parseETM5(strUTM, datum=_WGS84, Etm=Etm, falsed=True, **name):

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

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

1034

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

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

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

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

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

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

1041 @kwarg name: Optional B{C{Etm}} C{B{name}=NN} (C{str}).

1042

1043 @return: The UTM coordinate (B{C{Etm}}) or if C{B{Etm} is None}, a

1044 L{UtmUps5Tuple}C{(zone, hemipole, easting, northing, band)}

1045 with C{hemipole} is the hemisphere C{'N'|'S'}.

1046

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

1048

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

1050 '''

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

1052 return r

1053

1054

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

1056 strict=True, zone=None, **name_cmoff):

1057 '''Convert a geodetic lat-/longitude to an ETM coordinate.

1058

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

1060 C{LatLon} instance.

1061 @kwarg lon: Optional longitude (C{degrees}), required if B{C{latlon}}

1062 is C{degrees}, ignored otherwise.

1063 @kwarg datum: Optional datum for the ETM coordinate, overriding

1064 B{C{latlon}}'s datum (L{Datum}, L{Ellipsoid},

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

1066 @kwarg Etm: Optional class to return the ETM coordinate (L{Etm}) or C{None}.

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

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 name_cmoff: Optional B{C{Etm}} C{B{name}=NN} (C{str}) and DEPRECATED

1071 keyword argument C{B{cmoff}=True} to offset the longitude from

1072 the zone's central meridian (C{bool}), use B{C{falsed}} instead.

1073

1074 @return: The ETM coordinate as B{C{Etm}} or if C{B{Etm} is None} or not B{C{falsed}},

1075 a L{UtmUps8Tuple}C{(zone, hemipole, easting, northing, band, datum, gamma,

1076 scale)}. The C{hemipole} is the C{'N'|'S'} hemisphere.

1077

1078 @raise ETMError: No convergence transforming to ETM easting and northing.

1079

1080 @raise ETMError: Invalid B{C{zone}} or near-spherical or incompatible B{C{datum}}

1081 or C{ellipsoid}.

1082

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

1084 outside the valid range and L{rangerrors<pygeodesy.rangerrors>} is C{True}.

1085

1086 @raise TypeError: Invalid or near-spherical B{C{datum}} or B{C{latlon}} not ellipsoidal.

1087

1088 @raise ValueError: The B{C{lon}} value is missing or B{C{latlon}} is invalid.

1089 '''

1090 z, B, lat, lon, d, f, n = _to7zBlldfn(latlon, lon, datum,

1091 falsed, zone, strict,

1092 ETMError, **name_cmoff)

1093 lon0 = _cmlon(z) if f else None

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

1095

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

1097 n, latlon, d.exactTM, Error=ETMError)

1098

1099

1100 if __name__ == '__main__': # MCCABE 16

1101

1102 def _main():

1103

1104 from pygeodesy import fstr, KTransverseMercator

1105# from pygeodesy.interns import _DASH_ # from internals

1106 from pygeodesy.internals import printf, _usage, _DASH_

1107 from sys import argv, exit as _exit

1108

1109 def _help():

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

1111 '[-p[recision] <ndigits>',

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

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

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

1115 '|-h[elp]'))

1116

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

1118 _f = _r = _s = _t = False

1119 _p = -6

1120 _as = argv[1:]

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

1122 _a = _as.pop(0)

1123 if len(_a) < 2:

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

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

1126 _f, _r = True, False

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

1128 _f, _r = False, True

1129 elif '-precision'.startswith(_a):

1130 _p = int(_as.pop(0))

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

1132 _s, _t = True, False

1133 elif _a == '-t':

1134 _s, _t = False, True

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

1136 _help()

1137 else:

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

1139

1140 if len(_as) > 1:

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

1142 else:

1143 printf('%s ...: incomplete', _usage(*argv))

1144 _help()

1145

1146 if _s: # -series

1147 tm = KTransverseMercator()

1148 else:

1149 tm = ExactTransverseMercator(extendp=_t)

1150

1151 if _f:

1152 t = tm.forward(*f2)

1153 elif _r:

1154 t = tm.reverse(*f2)

1155 else:

1156 t = tm.forward(*f2)

1157 printf('%s: %s', tm.classname, fstr(t, prec=_p, sep=_SPACE_))

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

1159 printf('%s: %s', tm.classname, fstr(t, prec=_p, sep=_SPACE_))

1160

1161 _main()

1162

1163# % python3.13 -m pygeodesy.etm -p 12 33.33 44.44

1164# ExactTransverseMercator: 4276926.114803905599 4727193.767015309073 28.375536563148 1.233325101778

1165# ExactTransverseMercator: 33.33 44.44 28.375536563148 1.233325101778

1166

1167# % python3.13 -m pygeodesy.etm -s -p 12 33.33 44.44

1168# KTransverseMercator: 4276926.114803904667 4727193.767015310004 28.375536563148 1.233325101778

1169# KTransverseMercator: 33.33 44.44 28.375536563148 1.233325101778

1170

1171# % python3.12 -m pygeodesy.etm -p 12 33.33 44.44

1172# ExactTransverseMercator: 4276926.11480390653 4727193.767015309073 28.375536563148 1.233325101778

1173# ExactTransverseMercator: 33.33 44.44 28.375536563148 1.233325101778

1174

1175# % python3.12 -m pygeodesy.etm -s -p 12 33.33 44.44

1176# KTransverseMercator: 4276926.114803904667 4727193.767015310004 28.375536563148 1.233325101778

1177# KTransverseMercator: 33.33 44.44 28.375536563148 1.233325101778

1178

1179# % python2 -m pygeodesy.etm -p 12 33.33 44.44

1180# ExactTransverseMercator: 4276926.11480390653 4727193.767015309073 28.375536563148 1.233325101778

1181# ExactTransverseMercator: 33.33 44.44 28.375536563148 1.233325101778

1182

1183# % python2 -m pygeodesy.etm -s -p 12 33.33 44.44

1184# KTransverseMercator: 4276926.114803904667 4727193.767015310004 28.375536563148 1.233325101778

1185# KTransverseMercator: 33.33 44.44 28.375536563148 1.233325101778

1186

1187# % echo 33.33 44.44 | .../bin/TransverseMercatorProj

1188# 4276926.114804 4727193.767015 28.375536563148 1.233325101778

1189

1190# **) MIT License

1191#

1192# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1193#

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

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

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

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

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

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

1200#

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

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

1203#

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

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

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

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

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

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

1210# OTHER DEALINGS IN THE SOFTWARE.