1

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

3

4 u'''World Geographic Reference System (WGRS) en-/decoding.

5

6 Classes L{Georef} and L{WGRSError} and several functions to encode,

7 decode and inspect WGRS references.

8

9 Transcoded from I{Charles Karney}'s C++ class U{Georef

10<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Georef.html>},

11 but with modified C{precision} and extended with C{height} and C{radius}. See

12 also U{World Geographic Reference System

13<https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.

14'''

15# from pygeodesy.basics import isstr # from .named

16 from pygeodesy.constants import _float, _off90, _0_001, _0_5, \

17 _1_0, _2_0, _60_0, _1000_0

18 from pygeodesy.dms import parse3llh # parseDMS2

19 from pygeodesy.errors import _ValueError, _xattr, _xkwds

20 from pygeodesy.interns import NN, _0to9_, _AtoZnoIO_, _COMMA_, \

21 _height_, _radius_, _SPACE_

22 from pygeodesy.lazily import _ALL_LAZY, _ALL_OTHER

23 from pygeodesy.named import isstr, nameof, Property_RO

24 from pygeodesy.namedTuples import LatLon2Tuple, LatLonPrec3Tuple

25# from pygeodesy.props import Property_RO # from .named

26 from pygeodesy.streprs import Fmt, _0wd

27 from pygeodesy.units import Height, Int, Lat, Lon, Precision_, \

28 Radius, Scalar_, Str, _xStrError

29 from pygeodesy.utily import ft2m, m2ft, m2NM

30

31 from math import floor

32

33__all__ = _ALL_LAZY.wgrs

34__version__ = '23.05.06'

35

36_Base = 10

37_BaseLen = 4

38_DegChar = _AtoZnoIO_.tillQ

39_Digits = _0to9_

40_Height_ = Height.__name__

41_INV_ = 'INV'

42_LatOrig = -90

43_LatTile = _AtoZnoIO_.tillM

44_LonOrig = -180

45_LonTile = _AtoZnoIO_

46_60B = 60000000000 # == 60_000_000_000 == 60 * pow(10, 9)

47_MaxPrec = 11

48_Radius_ = Radius.__name__

49_Tile = 15 # tile size in degrees

50

51_MaxLen = _BaseLen + 2 * _MaxPrec

52_MinLen = _BaseLen - 2

53

54_LatOrig_60B = _LatOrig * _60B

55_LonOrig_60B = _LonOrig * _60B

56

57_float_Tile = _float(_Tile)

58_LatOrig_Tile = _float(_LatOrig) / _Tile

59_LonOrig_Tile = _float(_LonOrig) / _Tile

60

61

62 def _divmod3(x, _Orig_60B):

63 '''(INTERNAL) Convert B{C{x}} to 3_tuple C{(tile, modulo, fraction)}/

64 '''

65 i = int(floor(x * _60B))

66 i, x = divmod(i - _Orig_60B, _60B)

67 xt, xd = divmod(i, _Tile)

68 return xt, xd, x

69

70

71 def _fllh3(lat, lon, height=None):

72 '''(INTERNAL) Convert lat, lon, height.

73 '''

74 # lat, lon = parseDMS2(lat, lon)

75 return (Lat(lat, Error=WGRSError),

76 Lon(lon, Error=WGRSError), height)

77

78

79 def _geostr2(georef):

80 '''(INTERNAL) Check a georef string.

81 '''

82 try:

83 n, geostr = len(georef), georef.upper()

84 p, o = divmod(n, 2)

85 if o or n < _MinLen or n > _MaxLen \

86 or geostr[:3] == _INV_ \

87 or not geostr.isalnum():

88 raise ValueError

89 return geostr, _Precision(p - 1)

90

91 except (AttributeError, TypeError, ValueError) as x:

92 raise WGRSError(Georef.__name__, georef, cause=x)

93

94

95 def _Precision(precision):

96 '''(INTERNAL) Return a L{Precision_} instance.

97 '''

98 return Precision_(precision, Error=WGRSError, low=0, high=_MaxPrec)

99

100

101 class WGRSError(_ValueError):

102 '''World Geographic Reference System (WGRS) encode, decode or other L{Georef} issue.

103 '''

104 pass

105

106

107 class Georef(Str):

108 '''Georef class, a named C{str}.

109 '''

110 # no str.__init__ in Python 3

111 def __new__(cls, cll, precision=3, name=NN):

112 '''New L{Georef} from an other L{Georef} instance or georef

113 C{str} or from a C{LatLon} instance or lat-/longitude C{str}.

114

115 @arg cll: Cell or location (L{Georef} or C{str}, C{LatLon}

116 or C{str}).

117 @kwarg precision: Optional, the desired georef resolution

118 and length (C{int} 0..11), see function

119 L{wgrs.encode} for more details.

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

121

122 @return: New L{Georef}.

123

124 @raise RangeError: Invalid B{C{cll}} lat- or longitude.

125

126 @raise TypeError: Invalid B{C{cll}}.

127

128 @raise WGRSError: INValid or non-alphanumeric B{C{cll}}.

129 '''

130 ll = p = None

131

132 if isinstance(cll, Georef):

133 g, p = _geostr2(str(cll))

134

135 elif isstr(cll):

136 if _COMMA_ in cll:

137 lat, lon, h = _fllh3(*parse3llh(cll, height=None))

138 g = encode(lat, lon, precision=precision, height=h) # PYCHOK false

139 ll = lat, lon # original lat, lon

140 else:

141 g = cll.upper()

142

143 else: # assume LatLon

144 try:

145 lat, lon, h = _fllh3(cll.lat, cll.lon)

146 except AttributeError:

147 raise _xStrError(Georef, cll=cll) # Error=WGRSError

148 h = _xattr(cll, height=h)

149 g = encode(lat, lon, precision=precision, height=h) # PYCHOK false

150 ll = lat, lon # original lat, lon

151

152 self = Str.__new__(cls, g, name=name or nameof(cll))

153 self._latlon = ll

154 self._precision = p

155 return self

156

157 @Property_RO

158 def decoded3(self):

159 '''Get this georef's attributes (L{LatLonPrec3Tuple}).

160 '''

161 lat, lon = self.latlon

162 return LatLonPrec3Tuple(lat, lon, self.precision, name=self.name)

163

164 @Property_RO

165 def decoded5(self):

166 '''Get this georef's attributes (L{LatLonPrec5Tuple}) with

167 height and radius set to C{None} if missing.

168 '''

169 return self.decoded3.to5Tuple(self.height, self.radius)

170

171 @Property_RO

172 def _decoded5(self):

173 '''(INTERNAL) Initial L{LatLonPrec5Tuple}.

174 '''

175 return decode5(self)

176

177 @Property_RO

178 def height(self):

179 '''Get this georef's height in C{meter} or C{None} if missing.

180 '''

181 return self._decoded5.height

182

183 @Property_RO

184 def latlon(self):

185 '''Get this georef's (center) lat- and longitude (L{LatLon2Tuple}).

186 '''

187 lat, lon = self._latlon or self._decoded5[:2]

188 return LatLon2Tuple(lat, lon, name=self.name)

189

190 @Property_RO

191 def latlonheight(self):

192 '''Get this georef's (center) lat-, longitude and height (L{LatLon3Tuple}),

193 with height set to C{0} if missing.

194 '''

195 return self.latlon.to3Tuple(self.height or 0)

196

197 @Property_RO

198 def precision(self):

199 '''Get this georef's precision (C{int}).

200 '''

201 p = self._precision

202 return self._decoded5.precision if p is None else p

203

204 @Property_RO

205 def radius(self):

206 '''Get this georef's radius in C{meter} or C{None} if missing.

207 '''

208 return self._decoded5.radius

209

210 def toLatLon(self, LatLon=None, height=None, **LatLon_kwds):

211 '''Return (the center of) this georef cell as an instance

212 of the supplied C{LatLon} class.

213

214 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.

215 @kwarg height: Optional height (C{meter}).

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

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

218

219 @return: This georef location (B{C{LatLon}}) or if B{C{LatLon}}

220 is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)}.

221

222 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}}.

223 '''

224 if LatLon is None:

225 r = self.latlonheight if height is None else \

226 self.latlon.to3Tuple(height)

227 else:

228 h = (self.height or 0) if height is None else height

229 r = LatLon(*self.latlon, height=h, **_xkwds(LatLon_kwds, name=self.name))

230 return r

231

232

233 def decode3(georef, center=True):

234 '''Decode a C{georef} to lat-, longitude and precision.

235

236 @arg georef: To be decoded (L{Georef} or C{str}).

237 @kwarg center: If C{True} the center, otherwise the south-west,

238 lower-left corner (C{bool}).

239

240 @return: A L{LatLonPrec3Tuple}C{(lat, lon, precision)}.

241

242 @raise WGRSError: Invalid B{C{georef}}, INValid, non-alphanumeric

243 or odd length B{C{georef}}.

244 '''

245 def _digit(ll, g, i, m):

246 d = _Digits.find(g[i])

247 if d < 0 or d >= m:

248 raise _Error(i)

249 return ll * m + d

250

251 def _Error(i):

252 return WGRSError(Fmt.SQUARE(georef=i), georef)

253

254 def _index(chars, g, i):

255 k = chars.find(g[i])

256 if k < 0:

257 raise _Error(i)

258 return k

259

260 g, precision = _geostr2(georef)

261 lon = _index(_LonTile, g, 0) + _LonOrig_Tile

262 lat = _index(_LatTile, g, 1) + _LatOrig_Tile

263

264 u = _1_0

265 if precision > 0:

266 lon = lon * _Tile + _index(_DegChar, g, 2)

267 lat = lat * _Tile + _index(_DegChar, g, 3)

268 m, p = 6, precision - 1

269 for i in range(_BaseLen, _BaseLen + p):

270 lon = _digit(lon, g, i, m)

271 lat = _digit(lat, g, i + p, m)

272 u *= m

273 m = _Base

274 u *= _Tile

275

276 if center:

277 lon = lon * _2_0 + _1_0

278 lat = lat * _2_0 + _1_0

279 u *= _2_0

280 u = _Tile / u

281 return LatLonPrec3Tuple(Lat(lat * u, Error=WGRSError),

282 Lon(lon * u, Error=WGRSError),

283 precision, name=nameof(georef))

284

285

286 def decode5(georef, center=True):

287 '''Decode a C{georef} to lat-, longitude, precision, height and radius.

288

289 @arg georef: To be decoded (L{Georef} or C{str}).

290 @kwarg center: If C{True} the center, otherwise the south-west,

291 lower-left corner (C{bool}).

292

293 @return: A L{LatLonPrec5Tuple}C{(lat, lon,

294 precision, height, radius)} where C{height} and/or

295 C{radius} are C{None} if missing.

296

297 @raise WGRSError: Invalid B{C{georef}}, INValid, non-alphanumeric

298 or odd length B{C{georef}}.

299 '''

300 def _h2m(kft, name):

301 return Height(ft2m(kft * _1000_0), name=name, Error=WGRSError)

302

303 def _r2m(NM, name):

304 return Radius(NM / m2NM(1), name=name, Error=WGRSError)

305

306 def _split2(g, name, _2m):

307 i = max(g.find(name[0]), g.rfind(name[0]))

308 if i > _BaseLen:

309 return g[:i], _2m(int(g[i+1:]), _SPACE_(georef, name))

310 else:

311 return g, None

312

313 g = Str(georef, Error=WGRSError)

314

315 g, h = _split2(g, _Height_, _h2m) # H is last

316 g, r = _split2(g, _Radius_, _r2m) # R before H

317

318 return decode3(g, center=center).to5Tuple(h, r)

319

320

321 def encode(lat, lon, precision=3, height=None, radius=None): # MCCABE 14

322 '''Encode a lat-/longitude as a C{georef} of the given precision.

323

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

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

326 @kwarg precision: Optional, the desired C{georef} resolution and length

327 (C{int} 0..11).

328 @kwarg height: Optional, height in C{meter}, see U{Designation of area

329 <https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.

330 @kwarg radius: Optional, radius in C{meter}, see U{Designation of area

331 <https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.

332

333 @return: The C{georef} (C{str}).

334

335 @raise RangeError: Invalid B{C{lat}} or B{C{lon}}.

336

337 @raise WGRSError: Invalid B{C{precision}}, B{C{height}} or B{C{radius}}.

338

339 @note: The B{C{precision}} value differs from U{Georef<https://

340 GeographicLib.SourceForge.io/html/classGeographicLib_1_1Georef.html>}.

341 The C{georef} length is M{2 * (precision + 1)} and the

342 C{georef} resolution is I{15°} for B{C{precision}} 0, I{1°}

343 for 1, I{1′} for 2, I{0.1′} for 3, I{0.01′} for 4, ...

344 M{10**(2 - precision)}.

345 '''

346 def _option(name, m, m2_, K):

347 f = Scalar_(m, name=name, Error=WGRSError)

348 return NN(name[0].upper(), int(m2_(f * K) + _0_5))

349

350 p = _Precision(precision)

351

352 lat, lon, _ = _fllh3(lat, lon)

353 lat = _off90(lat)

354

355 xt, xd, x = _divmod3(lon, _LonOrig_60B)

356 yt, yd, y = _divmod3(lat, _LatOrig_60B)

357

358 g = _LonTile[xt], _LatTile[yt]

359 if p > 0:

360 g += _DegChar[xd], _DegChar[yd]

361 p -= 1

362 if p > 0:

363 d = pow(_Base, _MaxPrec - p)

364 x = _0wd(p, x // d)

365 y = _0wd(p, y // d)

366 g += x, y

367

368 if radius is not None: # R before H

369 g += _option(_radius_, radius, m2NM, _1_0),

370 if height is not None: # H is last

371 g += _option(_height_, height, m2ft, _0_001),

372

373 return NN.join(g) # XXX Georef(''.join(g))

374

375

376 def precision(res):

377 '''Determine the L{Georef} precision to meet a required (geographic)

378 resolution.

379

380 @arg res: The required resolution (C{degrees}).

381

382 @return: The L{Georef} precision (C{int} 0..11).

383

384 @raise ValueError: Invalid B{C{res}}.

385

386 @see: Function L{wgrs.encode} for more C{precision} details.

387 '''

388 r = Scalar_(res=res)

389 for p in range(_MaxPrec):

390 if resolution(p) <= r:

391 return p

392 return _MaxPrec

393

394

395 def resolution(prec):

396 '''Determine the (geographic) resolution of a given L{Georef} precision.

397

398 @arg prec: The given precision (C{int}).

399

400 @return: The (geographic) resolution (C{degrees}).

401

402 @raise ValueError: Invalid B{C{prec}}.

403

404 @see: Function L{wgrs.encode} for more C{precision} details.

405 '''

406 p = Int(prec=prec, Error=WGRSError)

407 if p > 1:

408 r = _1_0 / (_60_0 * pow(_Base, min(p, _MaxPrec) - 1))

409 elif p < 1:

410 r = _float_Tile

411 else:

412 r = _1_0

413 return r

414

415

416__all__ += _ALL_OTHER(decode3, decode5, # functions

417 encode, precision, resolution)

418

419# **) MIT License

420#

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

422#

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

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

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

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

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

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

429#

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

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

432#

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

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

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

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

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

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

439# OTHER DEALINGS IN THE SOFTWARE.