Coverage for pygeodesy/wgrs.py: 97%
189 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-08-02 18:24 -0400
« prev ^ index » next coverage.py v7.6.0, created at 2024-08-02 18:24 -0400
2# -*- coding: utf-8 -*-
4u'''World Geographic Reference System (WGRS) en-/decoding, aka GEOREF.
6Class L{Georef} and several functions to encode, decode and inspect
7WGRS (or GEOREF) references.
9Transcoded from I{Charles Karney}'s C++ class U{Georef
10<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Georef.html>},
11but with modified C{precision} and extended with C{height} and C{radius}.
13@see: U{World Geographic Reference System
14 <https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.
15'''
16# from pygeodesy.basics import isstr # from .named
17from pygeodesy.constants import INT0, _float, _off90, _0_001, \
18 _0_5, _1_0, _2_0, _60_0, _1000_0
19from pygeodesy.dms import parse3llh
20from pygeodesy.errors import _ValueError, _xattr, _xStrError
21from pygeodesy.interns import NN, _0to9_, _AtoZnoIO_, _COMMA_, \
22 _height_, _radius_, _SPACE_
23from pygeodesy.lazily import _ALL_LAZY, _ALL_OTHER
24from pygeodesy.named import _name2__, nameof, isstr, Property_RO
25from pygeodesy.namedTuples import LatLon2Tuple, LatLonPrec3Tuple
26# from pygeodesy.props import Property_RO # from .named
27from pygeodesy.streprs import Fmt, _0wd
28from pygeodesy.units import Height, Int, Lat, Lon, Precision_, \
29 Radius, Scalar_, Str
30from pygeodesy.utily import ft2m, m2ft, m2NM
32from math import floor
34__all__ = _ALL_LAZY.wgrs
35__version__ = '24.08.02'
37_Base = 10
38_BaseLen = 4
39_DegChar = _AtoZnoIO_.tillQ
40_Digits = _0to9_
41_INV_ = 'INV' # INValid
42_LatOrig = -90
43_LatTile = _AtoZnoIO_.tillM
44_LonOrig = -180
45_LonTile = _AtoZnoIO_
46_60B = 60000000000 # == 60_000_000_000 == 60e9
47_MaxPrec = 11
48_Tile = 15 # tile size in degrees
50_MaxLen = _BaseLen + 2 * _MaxPrec
51_MinLen = _BaseLen - 2
53_LatOrig_60B = _LatOrig * _60B
54_LonOrig_60B = _LonOrig * _60B
56_float_Tile = _float(_Tile)
57_LatOrig_Tile = _float(_LatOrig) / _Tile
58_LonOrig_Tile = _float(_LonOrig) / _Tile
61def _divmod3(x, _Orig_60B):
62 '''(INTERNAL) Convert B{C{x}} to 3_tuple C{(tile, modulo, fraction)}/
63 '''
64 i = int(floor(x * _60B))
65 i, x = divmod(i - _Orig_60B, _60B)
66 xt, xd = divmod(i, _Tile)
67 return xt, xd, x
70def _2fll(lat, lon):
71 '''(INTERNAL) Convert lat, lon.
72 '''
73 # lat, lon = parseDMS2(lat, lon)
74 return (Lat(lat, Error=WGRSError),
75 Lon(lon, Error=WGRSError))
78def _2geostr2(georef):
79 '''(INTERNAL) Check a georef string.
80 '''
81 try:
82 n, geostr = len(georef), georef.upper()
83 p, o = divmod(n, 2)
84 if o or n < _MinLen or n > _MaxLen \
85 or geostr[:3] == _INV_ \
86 or not geostr.isalnum():
87 raise ValueError()
88 return geostr, _2Precision(p - 1)
90 except (AttributeError, TypeError, ValueError) as x:
91 raise WGRSError(Georef.__name__, georef, cause=x)
94def _2Precision(precision):
95 '''(INTERNAL) Return a L{Precision_} instance.
96 '''
97 return Precision_(precision, Error=WGRSError, low=0, high=_MaxPrec)
100class WGRSError(_ValueError):
101 '''World Geographic Reference System (WGRS) encode, decode or other L{Georef} issue.
102 '''
103 pass
106class Georef(Str):
107 '''Georef class, a named C{str}.
108 '''
109 # no str.__init__ in Python 3
110 def __new__(cls, lat_gll, lon=None, height=None, precision=3, name=NN):
111 '''New L{Georef} from an other L{Georef} instance or georef
112 C{str} or from a C{LatLon} instance or lat-/longitude C{str}.
114 @arg lat_gll: Latitude (C{degrees90}), a georef (L{Georef},
115 C{str}) or a location (C{LatLon}, C{LatLon*Tuple}).
116 @kwarg lon: Logitude (C{degrees180)}, required if B{C{lat_gll}}
117 is C{degrees90}, ignored otherwise.
118 @kwarg height: Optional height in C{meter}, used if B{C{lat_gll}}
119 is a location.
120 @kwarg precision: The desired georef resolution and length (C{int}
121 0..11), see L{encode<pygeodesy.wgrs.encode>}.
122 @kwarg name: Optional name (C{str}).
124 @return: New L{Georef}.
126 @raise RangeError: Invalid B{C{lat_gll}} or B{C{lon}}.
128 @raise TypeError: Invalid B{C{lat_gll}} or B{C{lon}}.
130 @raise WGRSError: INValid B{C{lat_gll}}.
131 '''
132 if lon is None:
133 if isinstance(lat_gll, Georef):
134 g, ll, p = str(lat_gll), lat_gll.latlon, lat_gll.precision
135 elif isstr(lat_gll):
136 if _COMMA_ in lat_gll or _SPACE_ in lat_gll:
137 lat, lon, h = parse3llh(lat_gll, height=height)
138 g, ll, p = _encode3(lat, lon, precision, h=h)
139 else:
140 g, ll = lat_gll.upper(), None
141 try:
142 _, p = _2geostr2(g) # validate
143 except WGRSError: # R00H00?
144 p = None # = decode5(g).precision?
145 else: # assume LatLon
146 try:
147 g, ll, p = _encode3(lat_gll.lat, lat_gll.lon, precision,
148 h=_xattr(lat_gll, height=height))
149 except AttributeError:
150 raise _xStrError(Georef, gll=lat_gll) # Error=WGRSError
151 else:
152 g, ll, p = _encode3(lat_gll, lon, precision, h=height)
154 self = Str.__new__(cls, g, name=name or nameof(lat_gll))
155 self._latlon = ll
156 self._precision = p
157 return self
159 @Property_RO
160 def decoded3(self):
161 '''Get this georef's attributes (L{LatLonPrec3Tuple}).
162 '''
163 lat, lon = self.latlon
164 return LatLonPrec3Tuple(lat, lon, self.precision, name=self.name)
166 @Property_RO
167 def decoded5(self):
168 '''Get this georef's attributes (L{LatLonPrec5Tuple}) with
169 height and radius set to C{None} if missing.
170 '''
171 return self.decoded3.to5Tuple(self.height, self.radius)
173 @Property_RO
174 def _decoded5(self):
175 '''(INTERNAL) Initial L{LatLonPrec5Tuple}.
176 '''
177 return decode5(self)
179 @Property_RO
180 def height(self):
181 '''Get this georef's height in C{meter} or C{None} if missing.
182 '''
183 return self._decoded5.height
185 @Property_RO
186 def latlon(self):
187 '''Get this georef's (center) lat- and longitude (L{LatLon2Tuple}).
188 '''
189 lat, lon = self._latlon or self._decoded5[:2]
190 return LatLon2Tuple(lat, lon, name=self.name)
192 @Property_RO
193 def latlonheight(self):
194 '''Get this georef's (center) lat-, longitude and height (L{LatLon3Tuple}),
195 with height set to C{INT0} if missing.
196 '''
197 return self.latlon.to3Tuple(self.height or INT0)
199 @Property_RO
200 def precision(self):
201 '''Get this georef's precision (C{int}).
202 '''
203 p = self._precision
204 return self._decoded5.precision if p is None else p
206 @Property_RO
207 def radius(self):
208 '''Get this georef's radius in C{meter} or C{None} if missing.
209 '''
210 return self._decoded5.radius
212 def toLatLon(self, LatLon=None, height=None, **name_LatLon_kwds):
213 '''Return (the center of) this georef cell as a C{LatLon}.
215 @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
216 @kwarg height: Optional height (C{meter}), overriding this height.
217 @kwarg name_LatLon_kwds: Optional C{B{name}=NN} (C{str}) and optional,
218 additional B{C{LatLon}} keyword arguments, ignored if
219 C{B{LatLon} is None}.
221 @return: This georef location (B{C{LatLon}}) or if C{B{LatLon} is None},
222 a L{LatLon3Tuple}C{(lat, lon, height)}.
224 @raise TypeError: Invalid B{C{LatLon}} or B{C{name_LatLon_kwds}}.
225 '''
226 n, kwds = _name2__(name_LatLon_kwds, _or_nameof=self)
227 h = (self.height or INT0) if height is None else height # _heigHt
228 r = self.latlon.to3Tuple(h) if LatLon is None else LatLon(
229 *self.latlon, height=h, **kwds)
230 return r.renamed(n) if n else r
233def decode3(georef, center=True):
234 '''Decode a C{georef} to lat-, longitude and precision.
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}).
240 @return: A L{LatLonPrec3Tuple}C{(lat, lon, precision)}.
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
251 def _Error(i):
252 return WGRSError(Fmt.SQUARE(georef=i), georef)
254 def _index(chars, g, i):
255 k = chars.find(g[i])
256 if k < 0:
257 raise _Error(i)
258 return k
260 g, precision = _2geostr2(georef)
261 lon = _index(_LonTile, g, 0) + _LonOrig_Tile
262 lat = _index(_LatTile, g, 1) + _LatOrig_Tile
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
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))
286def decode5(georef, center=True):
287 '''Decode a C{georef} to lat-, longitude, precision, height and radius.
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}).
293 @return: A L{LatLonPrec5Tuple}C{(lat, lon, precision, height, radius)}
294 where C{height} and/or C{radius} are C{None} if missing.
296 @raise WGRSError: Invalid B{C{georef}}.
297 '''
298 def _h2m(kft, g_n):
299 return Height(ft2m(kft * _1000_0), name=g_n, Error=WGRSError)
301 def _r2m(NM, g_n):
302 return Radius(NM / m2NM(1), name=g_n, Error=WGRSError)
304 def _split2(g, Unit, _2m):
305 n = Unit.__name__
306 i = max(g.find(n[0]), g.rfind(n[0]))
307 if i > _BaseLen:
308 return g[:i], _2m(int(g[i+1:]), _SPACE_(georef, n))
309 else:
310 return g, None
312 g = Str(georef, Error=WGRSError)
314 g, h = _split2(g, Height, _h2m) # H is last
315 g, r = _split2(g, Radius, _r2m) # R before H
317 return decode3(g, center=center).to5Tuple(h, r)
320def encode(lat, lon, precision=3, height=None, radius=None): # MCCABE 14
321 '''Encode a lat-/longitude as a C{georef} of the given precision.
323 @arg lat: Latitude (C{degrees}).
324 @arg lon: Longitude (C{degrees}).
325 @kwarg precision: Optional, the desired C{georef} resolution and length
326 (C{int} 0..11).
327 @kwarg height: Optional, height in C{meter}, see U{Designation of area
328 <https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.
329 @kwarg radius: Optional, radius in C{meter}, see U{Designation of area
330 <https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.
332 @return: The C{georef} (C{str}).
334 @raise RangeError: Invalid B{C{lat}} or B{C{lon}}.
336 @raise WGRSError: Invalid B{C{precision}}, B{C{height}} or B{C{radius}}.
338 @note: The B{C{precision}} value differs from U{Georef<https://
339 GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Georef.html>}.
340 The C{georef} length is M{2 * (precision + 1)} and the
341 C{georef} resolution is I{15°} for B{C{precision}} 0, I{1°}
342 for 1, I{1′} for 2, I{0.1′} for 3, I{0.01′} for 4, ...
343 M{10**(2 - precision)}.
344 '''
345 def _option(name, m, m2_, K):
346 f = Scalar_(m, name=name, Error=WGRSError)
347 return NN(name[0].upper(), int(m2_(f * K) + _0_5))
349 g, _, _ = _encode3(lat, lon, precision)
350 if radius is not None: # R before H
351 g += _option(_radius_, radius, m2NM, _1_0)
352 if height is not None: # H is last
353 g += _option(_height_, height, m2ft, _0_001)
354 return g
357def _encode3(lat, lon, precision, h=None):
358 '''Return 3-tuple C{(georef, (lat, lon), p)}.
359 '''
360 p = _2Precision(precision)
362 lat, lon = _2fll(lat, lon)
364 if h is None:
365 xt, xd, x = _divmod3( lon, _LonOrig_60B)
366 yt, yd, y = _divmod3(_off90(lat), _LatOrig_60B)
368 g = _LonTile[xt], _LatTile[yt]
369 if p > 0:
370 g += _DegChar[xd], _DegChar[yd]
371 p -= 1
372 if p > 0:
373 d = pow(_Base, _MaxPrec - p)
374 x = _0wd(p, x // d)
375 y = _0wd(p, y // d)
376 g += x, y
377 g = NN.join(g)
378 else:
379 g = encode(lat, lon, precision=p, height=h)
381 return g, (lat, lon), p # XXX Georef(''.join(g))
384def precision(res):
385 '''Determine the L{Georef} precision to meet a required (geographic)
386 resolution.
388 @arg res: The required resolution (C{degrees}).
390 @return: The L{Georef} precision (C{int} 0..11).
392 @raise ValueError: Invalid B{C{res}}.
394 @see: Function L{wgrs.encode} for more C{precision} details.
395 '''
396 r = Scalar_(res=res)
397 for p in range(_MaxPrec):
398 if resolution(p) <= r:
399 return p
400 return _MaxPrec
403def resolution(prec):
404 '''Determine the (geographic) resolution of a given L{Georef} precision.
406 @arg prec: The given precision (C{int}).
408 @return: The (geographic) resolution (C{degrees}).
410 @raise ValueError: Invalid B{C{prec}}.
412 @see: Function L{wgrs.encode} for more C{precision} details.
413 '''
414 p = Int(prec=prec, Error=WGRSError)
415 if p > 1:
416 r = _1_0 / (_60_0 * pow(_Base, min(p, _MaxPrec) - 1))
417 elif p < 1:
418 r = _float_Tile
419 else:
420 r = _1_0
421 return r
424__all__ += _ALL_OTHER(decode3, decode5, # functions
425 encode, precision, resolution)
427# **) MIT License
428#
429# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
430#
431# Permission is hereby granted, free of charge, to any person obtaining a
432# copy of this software and associated documentation files (the "Software"),
433# to deal in the Software without restriction, including without limitation
434# the rights to use, copy, modify, merge, publish, distribute, sublicense,
435# and/or sell copies of the Software, and to permit persons to whom the
436# Software is furnished to do so, subject to the following conditions:
437#
438# The above copyright notice and this permission notice shall be included
439# in all copies or substantial portions of the Software.
440#
441# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
442# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
443# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
444# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
445# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
446# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
447# OTHER DEALINGS IN THE SOFTWARE.