1

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

3

4 u'''I{Veness}' Universal Transverse Mercator (UTM) projection.

5

6 Classes L{Utm} and L{UTMError} and functions L{parseUTM5}, L{toUtm8} and

7 L{utmZoneBand5}.

8

9 Pure Python implementation of UTM / WGS-84 conversion functions using

10 an ellipsoidal earth model, transcoded from JavaScript originals by

11 I{(C) Chris Veness 2011-2016} published under the same MIT Licence**, see

12 U{UTM<https://www.Movable-Type.co.UK/scripts/latlong-utm-mgrs.html>} and

13 U{Module utm<https://www.Movable-Type.co.UK/scripts/geodesy/docs/module-utm.html>}.

14

15 The U{UTM<https://WikiPedia.org/wiki/Universal_Transverse_Mercator_coordinate_system>}

16 system is a 2-dimensional Cartesian coordinate system providing another way

17 to identify locations on the surface of the earth. UTM is a set of 60

18 transverse Mercator projections, normally based on the WGS-84 ellipsoid.

19 Within each zone, coordinates are represented as B{C{easting}}s and B{C{northing}}s,

20 measured in metres.

21

22 This module includes some of I{Charles Karney}'s U{'Transverse Mercator with an

23 accuracy of a few nanometers'<https://ArXiv.org/pdf/1002.1417v3.pdf>}, 2011

24(building on Krüger's U{'Konforme Abbildung des Erdellipsoids in der Ebene'

25<https://bib.GFZ-Potsdam.DE/pub/digi/krueger2.pdf>}, 1912) and C++ class

26 U{TransverseMercator<https://GeographicLib.SourceForge.io/C++/doc/

27 classGeographicLib_1_1TransverseMercator.html>}.

28

29 Some other references are U{Universal Transverse Mercator coordinate system

30<https://WikiPedia.org/wiki/Universal_Transverse_Mercator_coordinate_system>},

31 U{Transverse Mercator Projection<https://GeographicLib.SourceForge.io/tm.html>}

32 and Henrik Seidel U{'Die Mathematik der Gauß-Krueger-Abbildung'

33<https://DE.WikiPedia.org/wiki/Gauß-Krüger-Koordinatensystem>}, 2006.

34'''

35

36 from pygeodesy.basics import len2, map2, neg # splice

37 from pygeodesy.constants import EPS, EPS0, _K0_UTM, _0_0, _0_0001

38 from pygeodesy.datums import _ellipsoidal_datum, _WGS84

39 from pygeodesy.dms import degDMS, parseDMS2

40 from pygeodesy.errors import MGRSError, RangeError, _ValueError, \

41 _xkwds_get

42 from pygeodesy.fmath import fdot3, hypot, hypot1

43 from pygeodesy.interns import MISSING, NN, _by_, _COMMASPACE_, _N_, \

44 _NS_, _outside_, _range_, _S_, _scale0_, \

45 _SPACE_, _UTM_, _V_, _X_, _zone_, _under

46 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

47# from pygeodesy.named import _xnamed # from .utmupsBase

48 from pygeodesy.namedTuples import EasNor2Tuple, UtmUps5Tuple, \

49 UtmUps8Tuple, UtmUpsLatLon5Tuple

50 from pygeodesy.props import deprecated_method, property_doc_, \

51 Property_RO

52 from pygeodesy.streprs import Fmt, unstr

53 from pygeodesy.units import Band, Int, Lat, Lon, Meter, Zone

54 from pygeodesy.utily import degrees90, degrees180, sincos2

55 from pygeodesy.utmupsBase import _hemi, _LLEB, _parseUTMUPS5, _to4lldn, \

56 _to3zBhp, _to3zll, _UPS_LATS, _UPS_ZONE, \

57 _UTM_LAT_MAX, _UTM_ZONE_MAX, \

58 _UTM_LAT_MIN, _UTM_ZONE_MIN, \

59 _UTM_ZONE_OFF_MAX, UtmUpsBase, _xnamed

60

61 from math import asinh, atan, atanh, atan2, cos, cosh, degrees, \

62 fabs, radians, sin, sinh, tan, tanh

63 from operator import mul as _mul

64

65__all__ = _ALL_LAZY.utm

66__version__ = '23.08.05'

67

68_Bands = 'CDEFGHJKLMNPQRSTUVWXX' # UTM latitude bands C..X (no

69# I|O) 8° each, covering 80°S to 84°N and X repeated for 80-84°N

70_bandLat_ = 'bandLat'

71_FalseEasting = Meter( 500e3) # falsed offset (C{meter})

72_FalseNorthing = Meter(10000e3) # falsed offset (C{meter})

73_SvalbardXzone = {32: 9, 34: 21, 36: 33} # [zone] longitude

74

75

76 class UTMError(_ValueError):

77 '''Universal Transverse Mercator (UTM parse or other L{Utm} issue.

78 '''

79 pass

80

81

82 class _Kseries(object):

83 '''(INTERNAL) Alpha or Beta Krüger series.

84

85 Krüger series summations for B{C{eta}}, B{C{ksi}}, B{C{p}} and B{C{q}},

86 caching the C{cos}, C{cosh}, C{sin} and C{sinh} values for

87 the given B{C{eta}} and B{C{ksi}} angles (in C{radians}).

88 '''

89 def __init__(self, AB, x, y):

90 '''(INTERNAL) New Alpha or Beta Krüger series

91

92 @arg AB: Krüger Alpha or Beta series coefficients

93 (C{4-, 6- or 8-tuple}).

94 @arg x: Eta angle (C{radians}).

95 @arg y: Ksi angle (C{radians}).

96 '''

97 n, j2 = len2(range(2, len(AB) * 2 + 1, 2))

98

99 self._ab = AB

100 self._pq = map2(_mul, j2, AB)

101# assert len(self._ab) == len(self._pq) == n

102

103 x2 = map2(_mul, j2, (x,) * n)

104 self._chx = map2(cosh, x2)

105 self._shx = map2(sinh, x2)

106# assert len(x2) == len(self._chx) == len(self._shx) == n

107

108 y2 = map2(_mul, j2, (y,) * n)

109 self._cy = map2(cos, y2)

110 self._sy = map2(sin, y2)

111 # self._sy, self._cy = splice(sincos2(*y2)) # PYCHOK false

112# assert len(y2) == len(self._cy) == len(self._sy) == n

113

114 def xs(self, x0):

115 '''(INTERNAL) Eta summation (C{float}).

116 '''

117 return fdot3(self._ab, self._cy, self._shx, start=x0)

118

119 def ys(self, y0):

120 '''(INTERNAL) Ksi summation (C{float}).

121 '''

122 return fdot3(self._ab, self._sy, self._chx, start=y0)

123

124 def ps(self, p0):

125 '''(INTERNAL) P summation (C{float}).

126 '''

127 return fdot3(self._pq, self._cy, self._chx, start=p0)

128

129 def qs(self, q0):

130 '''(INTERNAL) Q summation (C{float}).

131 '''

132 return fdot3(self._pq, self._sy, self._shx, start=q0)

133

134

135 def _cmlon(zone):

136 '''(INTERNAL) Central meridian longitude (C{degrees180}).

137 '''

138 return (zone * 6) - 183

139

140

141 def _false2(e, n, h):

142 '''(INTERNAL) False easting and northing.

143 '''

144 # Karney, "Test data for the transverse Mercator projection (2009)"

145 # <https://GeographicLib.SourceForge.io/C++/doc/transversemercator.html>

146 # and <https://Zenodo.org/record/32470#.W4LEJS2ZON8>

147 e += _FalseEasting # make e relative to central meridian

148 if h == _S_:

149 n += _FalseNorthing # make n relative to equator

150 return e, n

151

152

153 def _toBand(lat, *unused, **strict_Error): # see ups._toBand

154 '''(INTERNAL) Get the I{latitudinal} Band (row) letter.

155 '''

156 if _UTM_LAT_MIN <= lat < _UTM_LAT_MAX: # [-80, 84) like Veness

157 return _Bands[int(lat - _UTM_LAT_MIN) >> 3]

158 elif _xkwds_get(strict_Error, strict=True):

159 r = _range_(_UTM_LAT_MIN, _UTM_LAT_MAX, ropen=True)

160 t = _SPACE_(_outside_, _UTM_, _range_, r)

161 E = _xkwds_get(strict_Error, Error=RangeError)

162 raise E(lat=degDMS(lat), txt=t)

163 else:

164 return NN # None

165

166

167 def _to3zBlat(zone, band, Error=UTMError): # in .mgrs

168 '''(INTERNAL) Check and return zone, Band and band latitude.

169

170 @arg zone: Zone number or string.

171 @arg band: Band letter.

172 @arg Error: Exception to raise (L{UTMError}).

173

174 @return: 3-Tuple (zone, Band, latitude).

175 '''

176 z, B, _ = _to3zBhp(zone, band, Error=Error)

177 if not (_UTM_ZONE_MIN <= z <= _UTM_ZONE_MAX or

178 (_UPS_ZONE == z and Error is MGRSError)):

179 raise Error(zone=zone)

180

181 b = None

182 if B:

183 if z == _UPS_ZONE: # polar

184 try:

185 b = Lat(_UPS_LATS[B], name=_bandLat_)

186 except KeyError:

187 raise Error(band=band or B, zone=z)

188 else: # UTM

189 b = _Bands.find(B)

190 if b < 0:

191 raise Error(band=band or B, zone=z)

192 b = Int((b << 3) - 80, name=_bandLat_)

193 B = Band(B)

194 elif Error is not UTMError:

195 raise Error(band=band, txt=MISSING)

196

197 return Zone(z), B, b

198

199

200 def _to4zBll(lat, lon, cmoff=True, strict=True, Error=RangeError):

201 '''(INTERNAL) Return zone, Band and lat- and (central) longitude in degrees.

202

203 @arg lat: Latitude (C{degrees}).

204 @arg lon: Longitude (C{degrees}).

205 @kwarg cmoff: Offset B{C{lon}} from zone's central meridian.

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

207 @kwarg Error: Error for out of UTM range B{C{lat}}s.

208

209 @return: 4-Tuple (zone, Band, lat, lon).

210 '''

211 z, lat, lon = _to3zll(lat, lon) # in .utmupsBase

212

213 x = lon - _cmlon(z) # z before Norway/Svalbard

214 if fabs(x) > _UTM_ZONE_OFF_MAX:

215 t = _SPACE_(_outside_, _UTM_, _zone_, str(z), _by_, degDMS(x, prec=6))

216 raise Error(lon=degDMS(lon), txt=t)

217

218 B = _toBand(lat, strict=strict, Error=Error)

219 if B == _X_: # and 0 <= lon < 42: z = int(lon + 183) // 6 + 1

220 x = _SvalbardXzone.get(z, None)

221 if x: # Svalbard/Spitsbergen archipelago

222 z += 1 if lon >= x else -1

223 elif B == _V_ and z == 31 and lon >= 3:

224 z += 1 # SouthWestern Norway

225

226 if cmoff: # lon off central meridian

227 lon -= _cmlon(z) # z after Norway/Svalbard

228 return Zone(z), (Band(B) if B else None), Lat(lat), Lon(lon)

229

230

231 def _to7zBlldfn(latlon, lon, datum, falsed, name, zone, strict, Error, **cmoff):

232 '''(INTERNAL) Determine 7-tuple (zone, band, lat, lon, datum,

233 falsed, name) for methods L{toEtm8} and L{toUtm8}.

234 '''

235 f = falsed and _xkwds_get(cmoff, cmoff=True) # DEPRECATED

236 lat, lon, d, name = _to4lldn(latlon, lon, datum, name)

237 z, B, lat, lon = _to4zBll(lat, lon, cmoff=f, strict=strict)

238 if zone: # re-zone for ETM/UTM

239 r, _, _ = _to3zBhp(zone, B)

240 if r != z:

241 if not _UTM_ZONE_MIN <= r <= _UTM_ZONE_MAX:

242 raise Error(zone=zone)

243 if f: # re-offset from central meridian

244 lon += _cmlon(z) - _cmlon(r)

245 z = r

246 return z, B, lat, lon, d, f, name

247

248

249 class Utm(UtmUpsBase):

250 '''Universal Transverse Mercator (UTM) coordinate.

251 '''

252# _band = NN # latitudinal band letter ('C'|..|'X', no 'I'|'O')

253 _Bands = _Bands # latitudinal Band letters (C{tuple})

254 _Error = UTMError # or etm.ETMError

255# _scale = None # grid scale factor (C{scalar}) or C{None}

256 _scale0 = _K0_UTM # central scale factor (C{scalar})

257 _zone = 0 # longitudinal zone (C{int} 1..60)

258

259 def __init__(self, zone=31, hemisphere=_N_, easting=166022, # PYCHOK expected

260 northing=0, band=NN, datum=_WGS84, falsed=True,

261 gamma=None, scale=None, name=NN, **convergence):

262 '''New L{Utm} UTM coordinate.

263

264 @kwarg zone: Longitudinal UTM zone (C{int}, 1..60) or zone with/-out

265 I{latitudinal} Band letter (C{str}, '1C'|..|'60X').

266 @kwarg hemisphere: Northern or southern hemisphere (C{str}, C{'N[orth]'}

267 or C{'S[outh]'}).

268 @kwarg easting: Easting, see B{C{falsed}} (C{meter}).

269 @kwarg northing: Northing, see B{C{falsed}} (C{meter}).

270 @kwarg band: Optional, I{latitudinal} band (C{str}, 'C'|..|'X', no 'I'|'O').

271 @kwarg datum: Optional, this coordinate's datum (L{Datum}, L{Ellipsoid},

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

273 @kwarg falsed: If C{True}, both B{C{easting}} and B{C{northing}} are

274 falsed (C{bool}).

275 @kwarg gamma: Optional meridian convergence, bearing off grid North,

276 clockwise from true North (C{degrees}) or C{None}.

277 @kwarg scale: Optional grid scale factor (C{scalar}) or C{None}.

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

279 @kwarg convergence: DEPRECATED, use keyword argument C{B{gamma}=None}.

280

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

282

283 @raise UTMError: Invalid B{C{zone}}, B{C{hemishere}}, B{C{easting}},

284 B{C{northing}}, B{C{band}}, B{C{convergence}} or

285 B{C{scale}}.

286

287 @example:

288

289 >>> import pygeodesy

290 >>> u = pygeodesy.Utm(31, 'N', 448251, 5411932)

291 '''

292 if name:

293 self.name = name

294

295 self._zone, B, _ = _to3zBlat(zone, band)

296

297 h = str(hemisphere)[:1].upper()

298 if h not in _NS_:

299 raise self._Error(hemisphere=hemisphere)

300

301 e, n = easting, northing # Easting(easting), ...

302# if not falsed:

303# e, n = _false2(e, n, h)

304# # check easting/northing (with 40km overlap

305# # between zones) - is this worthwhile?

306# @raise RangeError: If B{C{easting}} or B{C{northing}} outside

307# the valid UTM range.

308# if 120e3 > e or e > 880e3:

309# raise RangeError(easting=easting)

310# if 0 > n or n > _FalseNorthing:

311# raise RangeError(northing=northing)

312

313 self._hemisphere = h

314 UtmUpsBase.__init__(self, e, n, band=B, datum=datum, falsed=falsed,

315 gamma=gamma, scale=scale, **convergence)

316

317 def __eq__(self, other):

318 return isinstance(other, Utm) and other.zone == self.zone \

319 and other.hemisphere == self.hemisphere \

320 and other.easting == self.easting \

321 and other.northing == self.northing \

322 and other.band == self.band \

323 and other.datum == self.datum

324

325 def __repr__(self):

326 return self.toRepr(B=True)

327

328 def __str__(self):

329 return self.toStr()

330

331 def _xcopy2(self, Xtm, name=NN):

332 '''(INTERNAL) Make copy as an B{C{Xtm}} instance.

333

334 @arg Xtm: Class to return the copy (C{Xtm=Etm}, C{Xtm=Utm} or

335 C{self.classof}).

336 '''

337 return Xtm(self.zone, self.hemisphere, self.easting, self.northing,

338 band=self.band, datum=self.datum, falsed=self.falsed,

339 gamma=self.gamma, scale=self.scale, name=name or self.name)

340

341 @property_doc_(''' the I{latitudinal} band.''')

342 def band(self):

343 '''Get the I{latitudinal} band (C{'C'|..|'X'}).

344 '''

345 if not self._band:

346 self._toLLEB()

347 return self._band

348

349 @band.setter # PYCHOK setter!

350 def band(self, band):

351 '''Set or reset the I{latitudinal} band letter (C{'C'|..|'X'})

352 or C{None} or C{""} to reset.

353

354 @raise TypeError: Invalid B{C{band}}.

355

356 @raise ValueError: Invalid B{C{band}}.

357 '''

358 self._band1(band)

359

360 @Property_RO

361 def _etm(self):

362 '''(INTERNAL) Cache for method L{toEtm}.

363 '''

364 return self._xcopy2(_MODS.etm.Etm)

365

366 @Property_RO

367 def falsed2(self):

368 '''Get the easting and northing falsing (L{EasNor2Tuple}C{(easting, northing)}).

369 '''

370 e = n = 0

371 if self.falsed:

372 e = _FalseEasting # relative to central meridian

373 if self.hemisphere == _S_: # relative to equator

374 n = _FalseNorthing

375 return EasNor2Tuple(e, n)

376

377 def parse(self, strUTM, name=NN):

378 '''Parse a string to a similar L{Utm} instance.

379

380 @arg strUTM: The UTM coordinate (C{str}),

381 see function L{parseUTM5}.

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

383 overriding this name.

384

385 @return: The similar instance (L{Utm}).

386

387 @raise UTMError: Invalid B{C{strUTM}}.

388

389 @see: Functions L{pygeodesy.parseUPS5} and L{pygeodesy.parseUTMUPS5}.

390 '''

391 return parseUTM5(strUTM, datum=self.datum, Utm=self.classof,

392 name=name or self.name)

393

394 @deprecated_method

395 def parseUTM(self, strUTM): # PYCHOK no cover

396 '''DEPRECATED, use method L{Utm.parse}.'''

397 return self.parse(strUTM)

398

399 @Property_RO

400 def pole(self):

401 '''Get the top center of (stereographic) projection, C{""} always.

402 '''

403 return NN # N/A for UTM

404

405 def toEtm(self):

406 '''Copy this UTM to an ETM coordinate.

407

408 @return: The ETM coordinate (L{Etm}).

409 '''

410 return self._etm

411

412 def toLatLon(self, LatLon=None, eps=EPS, unfalse=True, **LatLon_kwds):

413 '''Convert this UTM coordinate to an (ellipsoidal) geodetic point.

414

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

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

417 @kwarg eps: Optional convergence limit, L{EPS} or above (C{float}).

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

419 if falsed (C{bool}).

420 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword

421 arguments, ignored if C{B{LatLon} is None}.

422

423 @return: This UTM as (B{C{LatLon}}) or if B{C{LatLon}} is

424 C{None}, as L{LatLonDatum5Tuple}C{(lat, lon, datum,

425 gamma, scale)}.

426

427 @raise TypeError: Invalid B{C{datum}} or B{C{LatLon}} is not ellipsoidal.

428

429 @raise UTMError: Invalid meridional radius or H-value.

430

431 @example:

432

433 >>> u = Utm(31, 'N', 448251.795, 5411932.678)

434 >>> from pygeodesy import ellipsoidalVincenty as eV

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

436 '''

437 if eps < EPS:

438 eps = EPS # less doesn't converge

439

440 if self._latlon and self._latlon._toLLEB_args == (unfalse, eps):

441 return self._latlon5(LatLon)

442 else:

443 self._toLLEB(unfalse=unfalse, eps=eps)

444 return self._latlon5(LatLon, **LatLon_kwds)

445

446 def _toLLEB(self, unfalse=True, eps=EPS): # PYCHOK signature

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

448 '''

449 x, y = self.eastingnorthing2(falsed=not unfalse)

450

451 E = self.datum.ellipsoid

452 # from Karney 2011 Eq 15-22, 36

453 A0 = self.scale0 * E.A

454 if A0 < EPS0:

455 raise self._Error(meridional=A0)

456 x = x / A0 # /= chokes PyChecker

457 y = y / A0

458 K = _Kseries(E.BetaKs, x, y) # Krüger series

459 x = neg(K.xs(-x)) # η' eta

460 y = neg(K.ys(-y)) # ξ' ksi

461

462 sy, cy = sincos2(y)

463 shx = sinh(x)

464 H = hypot(shx, cy)

465 if H < EPS0:

466 raise self._Error(H=H)

467

468 T = sy / H # τʹ == τ0

469 p = _0_0 # previous d

470 e = _0_0001 * eps

471 for T, i, d in E._es_tauf3(T, T): # max 5

472 # d may toggle on +/-1.12e-16 or +/-4.47e-16,

473 # see the references at C{Ellipsoid.es_tauf}

474 if fabs(d) < eps or fabs(d + p) < e:

475 break

476 p = d

477 else:

478 t = unstr(self.toLatLon, eps=eps, unfalse=unfalse)

479 raise self._Error(Fmt.no_convergence(d), txt=t)

480

481 a = atan(T) # phi, lat

482 b = atan2(shx, cy)

483 if unfalse and self.falsed:

484 b += radians(_cmlon(self.zone))

485 ll = _LLEB(degrees90(a), degrees180(b), datum=self.datum, name=self.name)

486

487 # gamma and scale: Karney 2011 Eq 26, 27 and 28

488 p = neg(K.ps(-1))

489 q = K.qs(0)

490 s = hypot(p, q) * E.a / A0

491 ll._gamma = degrees(atan(tan(y) * tanh(x)) + atan2(q, p))

492 ll._scale = (E.e2s(sin(a)) * H * hypot1(T) / s) if s else s # INF?

493 ll._iteration = i

494 self._latlon5args(ll, _toBand, unfalse, eps)

495

496 def toRepr(self, prec=0, fmt=Fmt.SQUARE, sep=_COMMASPACE_, B=False, cs=False, **unused): # PYCHOK expected

497 '''Return a string representation of this UTM coordinate.

498

499 Note that UTM coordinates are rounded, not truncated (unlike

500 MGRS grid references).

501

502 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).

503 @kwarg fmt: Enclosing backets format (C{str}).

504 @kwarg sep: Optional separator between name:value pairs (C{str}).

505 @kwarg B: Optionally, include latitudinal band (C{bool}).

506 @kwarg cs: Optionally, include meridian convergence and grid

507 scale factor (C{bool} or non-zero C{int} to specify

508 the precison like B{C{prec}}).

509

510 @return: This UTM as a string C{"[Z:09[band], H:N|S, E:meter,

511 N:meter]"} plus C{", C:degrees, S:float"} if B{C{cs}} is

512 C{True} (C{str}).

513 '''

514 return self._toRepr(fmt, B, cs, prec, sep)

515

516 def toStr(self, prec=0, sep=_SPACE_, B=False, cs=False): # PYCHOK expected

517 '''Return a string representation of this UTM coordinate.

518

519 To distinguish from MGRS grid zone designators, a space is

520 left between the zone and the hemisphere.

521

522 Note that UTM coordinates are rounded, not truncated (unlike

523 MGRS grid references).

524

525 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).

526 @kwarg sep: Optional separator to join (C{str}) or C{None}

527 to return an unjoined C{tuple} of C{str}s.

528 @kwarg B: Optionally, include latitudinal band (C{bool}).

529 @kwarg cs: Optionally, include meridian convergence and grid

530 scale factor (C{bool} or non-zero C{int} to specify

531 the precison like B{C{prec}}).

532

533 @return: This UTM as a string with C{zone[band], hemisphere,

534 easting, northing, [convergence, scale]} in

535 C{"00 N|S meter meter"} plus C{" degrees float"} if

536 B{C{cs}} is C{True} (C{str}).

537

538 @example:

539

540 >>> u = Utm(3, 'N', 448251, 5411932.0001)

541 >>> u.toStr(4) # 03 N 448251.0 5411932.0001

542 >>> u.toStr(sep=', ') # 03 N, 448251, 5411932

543 '''

544

545 return self._toStr(self.hemisphere, B, cs, prec, sep)

546

547 def toUps(self, pole=NN, eps=EPS, falsed=True, **unused):

548 '''Convert this UTM coordinate to a UPS coordinate.

549

550 @kwarg pole: Optional top/center of the UPS projection,

551 (C{str}, 'N[orth]'|'S[outh]').

552 @kwarg eps: Optional convergence limit, L{EPS} or above

553 (C{float}), see method L{Utm.toLatLon}.

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

555

556 @return: The UPS coordinate (L{Ups}).

557 '''

558 u = self._ups

559 if u is None or u.pole != (pole or u.pole) or falsed != bool(u.falsed):

560 ll = self.toLatLon(LatLon=_LLEB, eps=eps, unfalse=True)

561 ups = _MODS.ups

562 self._ups = u = ups.toUps8(ll, Ups=ups.Ups, falsed=falsed,

563 name=self.name, pole=pole)

564 return u

565

566 def toUtm(self, zone, eps=EPS, falsed=True, **unused):

567 '''Convert this UTM coordinate to a different zone.

568

569 @arg zone: New UTM zone (C{int}).

570 @kwarg eps: Optional convergence limit, L{EPS} or above

571 (C{float}), see method L{Utm.toLatLon}.

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

573

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

575 '''

576 if zone == self.zone and falsed == self.falsed:

577 return self.copy()

578 elif zone:

579 u = self._utm

580 if u is None or u.zone != zone or falsed != u.falsed:

581 ll = self.toLatLon(LatLon=_LLEB, eps=eps, unfalse=True)

582 self._utm = u = toUtm8(ll, Utm=self.classof, falsed=falsed,

583 name=self.name, zone=zone)

584 return u

585 raise self._Error(zone=zone)

586

587 @Property_RO

588 def zone(self):

589 '''Get the (longitudinal) zone (C{int}, 1..60).

590 '''

591 return self._zone

592

593

594 def _parseUTM5(strUTM, datum, Xtm, falsed, Error=UTMError, name=NN): # imported by .etm

595 '''(INTERNAL) Parse a string representing a UTM coordinate,

596 consisting of C{"zone[band] hemisphere easting northing"},

597 see L{pygeodesy.parseETM5} and L{parseUTM5}.

598 '''

599 z, h, e, n, B = _parseUTMUPS5(strUTM, None, Error=Error)

600 if _UTM_ZONE_MIN > z or z > _UTM_ZONE_MAX or (B and B not in _Bands):

601 raise Error(strUTM=strUTM, zone=z, band=B)

602

603 if Xtm is None:

604 r = UtmUps5Tuple(z, h, e, n, B, Error=Error, name=name)

605 else:

606 r = Xtm(z, h, e, n, band=B, datum=datum, falsed=falsed)

607 if name:

608 r = _xnamed(r, name, force=True)

609 return r

610

611

612 def parseUTM5(strUTM, datum=_WGS84, Utm=Utm, falsed=True, name=NN):

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

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

615

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

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

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

619 @kwarg Utm: Optional class to return the UTM coordinate

620 (L{Utm}) or C{None}.

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

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

623

624 @return: The UTM coordinate (B{C{Utm}}) or if B{C{Utm}}

625 is C{None}, a L{UtmUps5Tuple}C{(zone, hemipole,

626 easting, northing, band)}. The C{hemipole} is

627 the C{'N'|'S'} hemisphere.

628

629 @raise UTMError: Invalid B{C{strUTM}}.

630

631 @raise TypeError: Invalid B{C{datum}}.

632

633 @example:

634

635 >>> u = parseUTM5('31 N 448251 5411932')

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

637 >>> u = parseUTM5('31 N 448251.8 5411932.7')

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

639 '''

640 r = _parseUTM5(strUTM, datum, Utm, falsed, name=name)

641 return r

642

643

644 def toUtm8(latlon, lon=None, datum=None, Utm=Utm, falsed=True,

645 name=NN, strict=True,

646 zone=None, **cmoff):

647 '''Convert a lat-/longitude point to a UTM coordinate.

648

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

650 geodetic C{LatLon} point.

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

652 @kwarg datum: Optional datum for this UTM coordinate,

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

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

655 @kwarg Utm: Optional class to return the UTM coordinate

656 (L{Utm}) or C{None}.

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

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

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

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

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

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

663

664 @return: The UTM coordinate (B{C{Utm}}) or if B{C{Utm}} is

665 C{None} or not B{C{falsed}}, a L{UtmUps8Tuple}C{(zone,

666 hemipole, easting, northing, band, datum, gamma,

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

668

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

670 B{C{lat}} or B{C{lon}} outside the valid range

671 and L{pygeodesy.rangerrors} set to C{True}.

672

673 @raise TypeError: Invalid B{C{datum}} or B{C{latlon}} not ellipsoidal.

674

675 @raise UTMError: Invalid B{C{zone}}.

676

677 @raise ValueError: If B{C{lon}} value is missing or if

678 B{C{latlon}} is invalid.

679

680 @note: Implements Karney’s method, using 8-th order Krüger series,

681 giving results accurate to 5 nm (or better) for distances

682 up to 3,900 Km from the central meridian.

683

684 @example:

685

686 >>> p = LatLon(48.8582, 2.2945) # 31 N 448251.8 5411932.7

687 >>> u = toUtm(p) # 31 N 448252 5411933

688 >>> p = LatLon(13.4125, 103.8667) # 48 N 377302.4 1483034.8

689 >>> u = toUtm(p) # 48 N 377302 1483035

690 '''

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

692 falsed, name, zone,

693 strict, UTMError, **cmoff)

694 d = _ellipsoidal_datum(d, name=name)

695 E = d.ellipsoid

696

697 a, b = radians(lat), radians(lon)

698 # easting, northing: Karney 2011 Eq 7-14, 29, 35

699 sb, cb = sincos2(b)

700

701 T = tan(a)

702 T12 = hypot1(T)

703 S = sinh(E.e * atanh(E.e * T / T12))

704

705 T_ = T * hypot1(S) - S * T12

706 H = hypot(T_, cb)

707

708 y = atan2(T_, cb) # ξ' ksi

709 x = asinh(sb / H) # η' eta

710

711 A0 = E.A * getattr(Utm, _under(_scale0_), _K0_UTM) # Utm is class or None

712

713 K = _Kseries(E.AlphaKs, x, y) # Krüger series

714 y = K.ys(y) * A0 # ξ

715 x = K.xs(x) * A0 # η

716

717 # convergence: Karney 2011 Eq 23, 24

718 p_ = K.ps(1)

719 q_ = K.qs(0)

720 g = degrees(atan(T_ / hypot1(T_) * tan(b)) + atan2(q_, p_))

721 # scale: Karney 2011 Eq 25

722 k = E.e2s(sin(a)) * T12 / H * (A0 / E.a * hypot(p_, q_))

723

724 return _toXtm8(Utm, z, lat, x, y,

725 B, d, g, k, f, name, latlon, EPS)

726

727

728 def _toXtm8(Xtm, z, lat, x, y, B, d, g, k, f, # PYCHOK 13+ args

729 name, latlon, eps, Error=UTMError):

730 '''(INTERNAL) Helper for methods L{toEtm8} and L{toUtm8}.

731 '''

732 h = _hemi(lat)

733 if f:

734 x, y = _false2(x, y, h)

735 if Xtm is None: # DEPRECATED

736 r = UtmUps8Tuple(z, h, x, y, B, d, g, k, Error=Error, name=name)

737 else:

738 r = _xnamed(Xtm(z, h, x, y, band=B, datum=d, falsed=f,

739 gamma=g, scale=k), name)

740 if isinstance(latlon, _LLEB) and d is latlon.datum: # see ups.toUtm8

741 r._latlon5args(latlon, _toBand, f, eps) # XXX weakref(latlon)?

742 latlon._gamma = g

743 latlon._scale = k

744 elif not r._band:

745 r._band = _toBand(lat)

746 return r

747

748

749 def utmZoneBand5(lat, lon, cmoff=False, name=NN):

750 '''Return the UTM zone number, Band letter, hemisphere and

751 (clipped) lat- and longitude for a given location.

752

753 @arg lat: Latitude in degrees (C{scalar} or C{str}).

754 @arg lon: Longitude in degrees (C{scalar} or C{str}).

755 @kwarg cmoff: Offset longitude from the zone's central

756 meridian (C{bool}).

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

758

759 @return: A L{UtmUpsLatLon5Tuple}C{(zone, band, hemipole,

760 lat, lon)} where C{hemipole} is the C{'N'|'S'}

761 UTM hemisphere.

762

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

764 B{C{lat}} or B{C{lon}} outside the valid range

765 and L{pygeodesy.rangerrors} set to C{True}.

766

767 @raise ValueError: Invalid B{C{lat}} or B{C{lon}}.

768 '''

769 lat, lon = parseDMS2(lat, lon)

770 z, B, lat, lon = _to4zBll(lat, lon, cmoff=cmoff)

771 return UtmUpsLatLon5Tuple(z, B, _hemi(lat), lat, lon, name=name)

772

773# **) MIT License

774#

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

776#

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

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

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

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

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

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

783#

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

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

786#

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

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

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

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

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

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

793# OTHER DEALINGS IN THE SOFTWARE.