Coverage for pygeodesy/etm.py: 97%

401 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-11 14:35 -0400

1 

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

3 

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

5 

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

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

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

9 

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

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

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

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

14the results. 

15 

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

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

18file C{Header}. 

19 

20Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2022) and licensed 

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

22GeographicLib.SourceForge.io>} documentation. 

23 

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

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

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

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

28suitable starting guess. 

29 

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

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

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

33exceptions:: 

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 

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

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

61which 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 

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

65 

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

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

68 _1_EPS, isnear0, _0_0, _0_1, _0_5, \ 

69 _1_0, _2_0, _3_0, _4_0, _90_0, _180_0 

70from pygeodesy.datums import _ellipsoidal_datum, _WGS84 

71from pygeodesy.elliptic import _ALL_LAZY, Elliptic 

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

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

74from pygeodesy.fsums import Fsum, fsum1_ 

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

76 _spherical_, _usage 

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

78 _tand, _unsigned2 

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

80from pygeodesy.named import callername, _incompatible, _NamedBase 

81from pygeodesy.namedTuples import Forward4Tuple, Reverse4Tuple 

82from pygeodesy.props import deprecated_method, deprecated_property_RO, \ 

83 Property_RO, property_RO, _update_all, \ 

84 property_doc_ 

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

86from pygeodesy.units import Degrees, Scalar_ 

87from pygeodesy.utily import atand, atan2d, sincos2 

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

89 _to7zBlldfn, Utm, UTMError 

90 

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

92 

93__all__ = _ALL_LAZY.etm 

94__version__ = '22.10.12' 

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 

103def _overflow(x): 

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

105 ''' 

106 return _copyBit(_OVERFLOW, x) 

107 

108 

109class ETMError(UTMError): 

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

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

112 ''' 

113 pass 

114 

115 

116class Etm(Utm): 

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

118 a Universal Transverse Mercator (UTM) coordinate using the 

119 L{ExactTransverseMercator} projection for highest accuracy. 

120 

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

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

123 

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

125 html/classGeographicLib_1_1TransverseMercatorExact.html#details>}. 

126 ''' 

127 _Error = ETMError # see utm.UTMError 

128 _exactTM = None 

129 

130 __init__ = Utm.__init__ 

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

132 

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

134 

135 @example: 

136 

137 >>> import pygeodesy 

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

139 ''' 

140 

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

142 def exactTM(self): 

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

144 ''' 

145 if self._exactTM is None: 

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

147 return self._exactTM 

148 

149 @exactTM.setter # PYCHOK setter! 

150 def exactTM(self, exactTM): 

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

152 

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

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

155 ''' 

156 _xinstanceof(ExactTransverseMercator, exactTM=exactTM) 

157 

158 E = self.datum.ellipsoid 

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

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

161 self._exactTM = exactTM 

162 self._scale0 = exactTM.k0 

163 

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

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

166 

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

168 see function L{parseETM5}. 

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

170 overriding this name. 

171 

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

173 

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

175 

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

177 and L{pygeodesy.parseUTMUPS5}. 

178 ''' 

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

180 name=name or self.name) 

181 

182 @deprecated_method 

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

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

185 ''' 

186 return self.parse(strETM) 

187 

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

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

190 

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

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

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

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

195 

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

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

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

199 

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

201 incompatible or no convergence transforming to 

202 lat- and longitude. 

203 

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

205 

206 @example: 

207 

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

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

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

211 ''' 

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

213 self._toLLEB(unfalse=unfalse) 

214 return self._latlon5(LatLon) 

215 

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

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

218 ''' 

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

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

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

222 t = repr(d.ellipsoid) 

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

224 

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

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

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

228 

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

230 ll._gamma = g 

231 ll._scale = k 

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

233 

234 def toUtm(self): # PYCHOK signature 

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

236 

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

238 ''' 

239 return self._xcopy2(Utm) 

240 

241 

242class ExactTransverseMercator(_NamedBase): 

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

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

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

246 ''' 

247 _datum = None # Datum 

248 _E = None # Ellipsoid 

249 _extendp = False # use extended domain 

250# _iteration = None # ._sigmaInv2 and ._zetaInv2 

251 _k0 = _K0_UTM # central scale factor 

252 _lon0 = _0_0 # central meridian 

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

254 _mv = _1_0 # _1_0 - ._mu 

255 _raiser = False # throw Error 

256 _sigmaC = None # _sigmaInv04 case 

257 _zetaC = None # _zetaInv04 case 

258 

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

260 '''New L{ExactTransverseMercator} projection. 

261 

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

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

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

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

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

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

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

269 

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

271 or B{C{k0}}. 

272 

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

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

275 especially on B{X{extendp}}. 

276 

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

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

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

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

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

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

283 ''' 

284 if extendp: 

285 self._extendp = True 

286 if name: 

287 self.name = name 

288 if raiser: 

289 self.raiser = True 

290 

291 self.datum = datum # invokes ._reset 

292 self.k0 = k0 

293 self.lon0 = lon0 

294 

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

296 def datum(self): 

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

298 ''' 

299 return self._datum 

300 

301 @datum.setter # PYCHOK setter! 

302 def datum(self, datum): 

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

304 

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

306 ''' 

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

308 self._reset(d) 

309 self._datum = d 

310 

311 @Property_RO 

312 def _e(self): 

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

314 ''' 

315 return self._E.e 

316 

317 @Property_RO 

318 def _1_e_90(self): # PYCHOK no cover 

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

320 ''' 

321 return (_1_0 - self._e) * _90_0 

322 

323 @property_RO 

324 def ellipsoid(self): 

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

326 ''' 

327 return self._E 

328 

329 @Property_RO 

330 def _e_PI_2(self): 

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

332 ''' 

333 return self._e * PI_2 

334 

335 @Property_RO 

336 def _e_PI_4_(self): 

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

338 ''' 

339 return -self._e * PI_4 

340 

341 @Property_RO 

342 def _1_e_PI_2(self): 

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

344 ''' 

345 return (_1_0 - self._e) * PI_2 

346 

347 @Property_RO 

348 def _1_2e_PI_2(self): 

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

350 ''' 

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

352 

353 @property_RO 

354 def equatoradius(self): 

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

356 ''' 

357 return self._E.a 

358 

359 a = equatoradius 

360 

361 @Property_RO 

362 def _e_TAYTOL(self): 

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

364 ''' 

365 return self._e * _TAYTOL 

366 

367 @Property_RO 

368 def _Eu(self): 

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

370 ''' 

371 return Elliptic(self._mu) 

372 

373 @Property_RO 

374 def _Eu_cE(self): 

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

376 ''' 

377 return self._Eu.cE 

378 

379 def _Eu_2cE_(self, xi): 

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

381 ''' 

382 return self._Eu_cE * _2_0 - xi 

383 

384 @Property_RO 

385 def _Eu_cE_4(self): 

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

387 ''' 

388 return self._Eu_cE / _4_0 

389 

390 @Property_RO 

391 def _Eu_cK(self): 

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

393 ''' 

394 return self._Eu.cK 

395 

396 @Property_RO 

397 def _Eu_cK_cE(self): 

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

399 ''' 

400 return self._Eu_cK / self._Eu_cE 

401 

402 @Property_RO 

403 def _Eu_2cK_PI(self): 

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

405 ''' 

406 return self._Eu_cK / PI_2 

407 

408 @Property_RO 

409 def _Ev(self): 

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

411 ''' 

412 return Elliptic(self._mv) 

413 

414 @Property_RO 

415 def _Ev_cK(self): 

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

417 ''' 

418 return self._Ev.cK 

419 

420 @Property_RO 

421 def _Ev_cKE(self): 

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

423 ''' 

424 return self._Ev.cKE 

425 

426 @Property_RO 

427 def _Ev_3cKE_4(self): 

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

429 ''' 

430 return self._Ev_cKE * 0.75 

431 

432 @Property_RO 

433 def _Ev_5cKE_4(self): 

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

435 ''' 

436 return self._Ev_cKE * 1.25 

437 

438 @Property_RO 

439 def extendp(self): 

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

441 ''' 

442 return self._extendp 

443 

444 @property_RO 

445 def flattening(self): 

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

447 ''' 

448 return self._E.f 

449 

450 f = flattening 

451 

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

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

454 

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

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

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

458 the default if not C{None}. 

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

460 

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

462 

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

464 real &x, real &y, 

465 real &gamma, real &k)}. 

466 

467 @raise ETMError: No convergence, thrown iff property 

468 C{B{raiser}=True}. 

469 ''' 

470 lat = _fix90(lat) 

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

472 if self.extendp: 

473 backside = _lat = _lon = False 

474 else: # enforce the parity 

475 lat, _lat = _unsigned2(lat) 

476 lon, _lon = _unsigned2(lon) 

477 backside = lon > 90 

478 if backside: # PYCHOK no cover 

479 lon = _180_0 - lon 

480 if lat == 0: 

481 _lat = True 

482 

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

484 if lat == 90: 

485 u = self._Eu_cK 

486 v = self._iteration = self._zetaC = 0 

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

488 u = self._iteration = self._zetaC = 0 

489 v = self._Ev_cK 

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

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

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

493 

494 sncndn6 = self._sncndn6(u, v) 

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

496 g, k = self._zetaScaled(sncndn6, ll=False) \ 

497 if lat != 90 else (lon, self.k0) 

498 

499 if backside: 

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

501 y *= self._k0_a 

502 x *= self._k0_a 

503 if _lat: 

504 y, g = neg_(y, g) 

505 if _lon: 

506 x, g = neg_(x, g) 

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

508 name=name or self.name) 

509 

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

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

512 ''' 

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

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

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

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

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

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

519 s, c = sincos2(a) 

520 h = hypot(psi, dlam) 

521 r = cbrt(h * _3_mv_e) 

522 u = r * c 

523 v = r * s + self._Ev_cK 

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

525 return u, v, h 

526 

527 @property_RO 

528 def iteration(self): 

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

530 or C{ExactTransverseMercator.reverse} iteration number 

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

532 ''' 

533 return self._iteration 

534 

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

536 def k0(self): 

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

538 ''' 

539 return self._k0 # aka scale0 

540 

541 @k0.setter # PYCHOK setter! 

542 def k0(self, k0): 

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

544 

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

546 ''' 

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

548 if self._k0 != k0: 

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

550 self._k0 = k0 

551 

552 @Property_RO 

553 def _k0_a(self): 

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

555 ''' 

556 return self.k0 * self.equatoradius 

557 

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

559 def lon0(self): 

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

561 ''' 

562 return self._lon0 

563 

564 @lon0.setter # PYCHOK setter! 

565 def lon0(self, lon0): 

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

567 

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

569 ''' 

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

571 

572 @deprecated_property_RO 

573 def majoradius(self): # PYCHOK no cover 

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

575 return self.equatoradius 

576 

577 @Property_RO 

578 def _1_mu_2(self): 

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

580 ''' 

581 return _1_0 + self._mu * _0_5 

582 

583 @Property_RO 

584 def _3_mv(self): 

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

586 ''' 

587 return _3_0 / self._mv 

588 

589 @Property_RO 

590 def _3_mv_e(self): 

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

592 ''' 

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

594 

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

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

597 

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

599 

600 @raise ETMError: No convergence. 

601 ''' 

602 sca1, tol2 = _1_0, _TOL_10 

603 if psi: # _zetaInv2 

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

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

606 

607 _zeta3 = self._zeta3 

608 _zetaDwd2 = self._zetaDwd2 

609 else: # _sigmaInv2 

610 _zeta3 = self._sigma3 

611 _zetaDwd2 = self._sigmaDwd2 

612 

613 d2, r = tol2, self.raiser 

614 _U_2_ = Fsum(u).fsum2_ 

615 _V_2_ = Fsum(v).fsum2_ 

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

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

618 sncndn6 = self._sncndn6(u, v) 

619 du, dv = _zetaDwd2(*sncndn6) 

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

621 T = (taup - T) * sca1 

622 L -= lam 

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

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

625 if d2 < tol2: 

626 r = False 

627 break 

628 d2 = hypot2(dU, dV) 

629 

630 self._iteration = i 

631 if r: # PYCHOK no cover 

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

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

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

635 return u, v 

636 

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

638 def raiser(self): 

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

640 ''' 

641 return self._raiser 

642 

643 @raiser.setter # PYCHOK setter! 

644 def raiser(self, raiser): 

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

646 for convergence failures. 

647 ''' 

648 self._raiser = bool(raiser) 

649 

650 def _reset(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, name=NN): 

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: Central meridian (C{degrees180}), overriding 

678 the default if not C{None}. 

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

680 

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

682 

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

684 real &lat, real &lon, 

685 real &gamma, real &k)} 

686 

687 @raise ETMError: No convergence, thrown iff property 

688 C{B{raiser}=True}. 

689 ''' 

690 # undoes the steps in .forward. 

691 xi = y / self._k0_a 

692 eta = x / self._k0_a 

693 if self.extendp: 

694 backside = _lat = _lon = False 

695 else: # enforce the parity 

696 eta, _lon = _unsigned2(eta) 

697 xi, _lat = _unsigned2(xi) 

698 backside = xi > self._Eu_cE 

699 if backside: # PYCHOK no cover 

700 xi = self._Eu_2cE_(xi) 

701 

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

703 if xi or eta != self._Ev_cKE: 

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

705 else: # PYCHOK no cover 

706 u = self._iteration = self._sigmaC = 0 

707 v = self._Ev_cK 

708 

709 if v or u != self._Eu_cK: 

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

711 else: # PYCHOK no cover 

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

713 

714 if backside: # PYCHOK no cover 

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

716 if _lat: 

717 lat, g = neg_(lat, g) 

718 if _lon: 

719 lon, g = neg_(lon, g) 

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

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

722 iteration=self._iteration, 

723 name=name or self.name) 

724 

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

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

727 

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

729 from C{._zeta3}. 

730 

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

732 

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

734 real snu, real cnu, real dnu, 

735 real snv, real cnv, real dnv, 

736 real &gamma, real &k)}. 

737 ''' 

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

739 cnudnv = cnu * dnv 

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

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

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

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

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

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

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

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

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

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

750 if d2 > 0: 

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

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

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

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

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

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

757 k2 = fsum1_(mu, mv, mv * tau**2) 

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

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

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

761 else: 

762 k = _OVERFLOW 

763 return g, k 

764 

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

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

767 

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

769 

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

771 real v, real snv, real cnv, real dnv, 

772 real &xi, real &eta)}. 

773 

774 @raise ETMError: No convergence. 

775 ''' 

776 mu = self._mu * cnu 

777 mv = self._mv * cnv 

778 # Lee 55.4 writing 

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

780 d2 = cnu * mu + cnv * mv 

781 mu *= snu * dnu 

782 mv *= snv * dnv 

783 if d2 > 0: # /= chokes PyChecker 

784 mu = mu / d2 

785 mv = mv / d2 

786 else: 

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

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

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

790 return xi, v, d2 

791 

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

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

794 

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

796 

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

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

799 real &du, real &dv)}. 

800 ''' 

801 snuv = snu * snv 

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

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

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

805 r = cnv * dnu * dnv 

806 i = cnu * snuv * self._mu 

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

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

809 return du, dv 

810 

811 def _sigmaInv2(self, xi, eta): 

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

813 

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

815 

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

817 real &u, real &v)}. 

818 

819 @raise ETMError: No convergence. 

820 ''' 

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

822 if not t: 

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

824 return u, v 

825 

826 def _sigmaInv04(self, xi, eta): 

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

828 

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

830 

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

832 real &u, real &v)}. 

833 ''' 

834 t = False 

835 d = eta - self._Ev_cKE 

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

837 # sigma as a simple pole at 

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

839 # and sigma is approximated by 

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

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

842 u += self._Eu_cK 

843 v += self._Ev_cK 

844 C = 1 

845 

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

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

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

849 # sigma' = sigma'' = 0 

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

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

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

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

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

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

856 t = h < _TAYTOL2 

857 C = 2 

858 

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

860 u = v = self._Eu_cK_cE 

861 u *= xi 

862 v *= eta 

863 C = 3 

864 

865 return u, v, t, C 

866 

867 def _sncndn6(self, u, v): 

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

869 ''' 

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

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

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

873 

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

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

876 

877 @kwarg joined: Separator to join the attribute strings 

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

879 @kwarg kwds: Optional, overriding keyword arguments. 

880 ''' 

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

882 k0=self.k0, extendp=self.extendp) 

883 if self.name: 

884 d.update(name=self.name) 

885 t = pairs(d, **kwds) 

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

887 

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

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

890 

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

892 

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

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

895 real &taup, real &lam)} 

896 ''' 

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

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

899 t1 = t2 = _overflow(snu) 

900 # Lee 54.17 but write 

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

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

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

904 if d1 > EPS02: # _EPSmin 

905 t1 = snu * dnv / sqrt(d1) 

906 else: 

907 d1 = 0 

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

909 if d2 > EPS02: # _EPSmin 

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

911 else: 

912 d2 = 0 

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

914 # taup = sinh(psi) 

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

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

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

918 return taup, lam, d2 

919 

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

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

922 

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

924 

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

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

927 real &du, real &dv)}. 

928 ''' 

929 cnu2 = cnu**2 * self._mu 

930 cnv2 = cnv**2 

931 dnuv = dnu * dnv 

932 dnuv2 = dnuv**2 

933 snuv = snu * snv 

934 snuv2 = snuv**2 * self._mu 

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

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

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

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

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

940 return du, neg(dv) 

941 

942 def _zetaInv2(self, taup, lam): 

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

944 

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

946 

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

948 real &u, real &v)}. 

949 

950 @raise ETMError: No convergence. 

951 ''' 

952 psi = asinh(taup) 

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

954 if not t: 

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

956 return u, v 

957 

958 def _zetaInv04(self, psi, lam): 

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

960 

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

962 

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

964 real &u, real &v)}. 

965 ''' 

966 if lam > self._1_2e_PI_2: 

967 d = lam - self._1_e_PI_2 

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

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

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

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

972 # pole, where we have, approximately 

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

974 # Inverting this gives: 

975 e = self._e # eccentricity 

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

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

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

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

980 return u, v, False, 1 

981 

982 elif psi < self._e_PI_2: 

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

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

985 # zeta' = zeta'' = 0 

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

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

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

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

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

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

992 

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

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

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

996 # corresponding to the north pole. 

997 s, c = sincos2(lam) 

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

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

1000 u = r * atan2(h, c) 

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

1002 return u, v, False, 3 

1003 

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

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

1006 

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

1008 

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

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

1011 ''' 

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

1013 tau = self._E.es_tauf(t) 

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

1015 if ll: 

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

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

1018 

1019 

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

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

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

1023 

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

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

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

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

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

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

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

1031 

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

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

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

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

1036 

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

1038 

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

1040 

1041 @example: 

1042 

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

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

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

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

1047 ''' 

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

1049 return r 

1050 

1051 

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

1053 name=NN, strict=True, 

1054 zone=None, **cmoff): 

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

1056 

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

1058 geodetic C{LatLon} point. 

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

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

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

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

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

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

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

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

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

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

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

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

1071 

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

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

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

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

1076 hemisphere. 

1077 

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

1079 and northing. 

1080 

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

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

1083 

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

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

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

1087 

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

1089 B{C{latlon}} not ellipsoidal. 

1090 

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

1092 is invalid. 

1093 ''' 

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

1095 falsed, name, zone, 

1096 strict, ETMError, **cmoff) 

1097 lon0 = _cmlon(z) if f else None 

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

1099 

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

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

1102 

1103 

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

1105 

1106 from pygeodesy.lazily import _ALL_MODS as _MODS, printf 

1107 from sys import argv, exit as _exit 

1108 

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

1110 _f = _r = _s = _t = False 

1111 _as = argv[1:] 

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

1113 _a = _as.pop(0) 

1114 if len(_a) < 2: 

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

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

1117 _f, _r = True, False 

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

1119 _f, _r = False, True 

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

1121 _s, _t = True, False 

1122 elif _a == '-t': 

1123 _s, _t = False, True 

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

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

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

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

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

1129 '| -h[elp]')) 

1130 else: 

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

1132 if len(_as) > 1: 

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

1134 else: 

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

1136 

1137 if _s: 

1138 tm = _MODS.ktm.KTransverseMercator() 

1139 else: 

1140 tm = ExactTransverseMercator(extendp=_t) 

1141 

1142 if _f: 

1143 t = tm.forward(*f2) 

1144 elif _r: 

1145 t = tm.reverse(*f2) 

1146 else: 

1147 t = tm.forward(*f2) 

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

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

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

1151 

1152# **) MIT License 

1153# 

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

1155# 

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

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

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

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

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

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

1162# 

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

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

1165# 

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

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

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

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

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

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

1172# OTHER DEALINGS IN THE SOFTWARE.