1
2# -*- coding: utf-8 -*-
3
4u'''World Geographic Reference System (WGRS) en-/decoding.
5
6Classes L{Georef} and L{WGRSError} and several functions to encode,
7decode and inspect WGRS references.
8
9Transcoded from I{Charles Karney}'s C++ class U{Georef
10<https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1Georef.html>},
11but with modified C{precision} and extended with C{height} and C{radius}.  See
12also U{World Geographic Reference System
13<https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.
14'''
15
16from pygeodesy.basics import isstr
17from pygeodesy.dms import parse3llh  # parseDMS2
18from pygeodesy.errors import _ValueError, _xkwds
19from pygeodesy.interns import EPS1_2, NN, _AtoZnoIO_, _float, \
20                             _height_, _radius_, _SPACE_, \
21                             _0to9_, _0_5, _0_001, _1_0, _2_0, \
22                             _60_0, _90_0, _1000_0
23from pygeodesy.lazily import _ALL_LAZY, _ALL_OTHER
24from pygeodesy.named import nameof
25from pygeodesy.namedTuples import LatLon2Tuple, LatLonPrec3Tuple
26from pygeodesy.props import Property_RO
27from pygeodesy.streprs import Fmt, _0wd
28from pygeodesy.units import Height, Int, Lat, Lon, Precision_, \
29                            Radius, Scalar_, Str, _xStrError
30from pygeodesy.utily import ft2m, m2ft, m2NM
31
32from math import floor
33
34__all__ = _ALL_LAZY.wgrs
35__version__ = '21.07.31'
36
37_Base    =  10
38_BaseLen =  4
39_DegChar = _AtoZnoIO_.tillQ
40_Digits  = _0to9_
41_Height_ =  Height.__name__
42_INV_    = 'INV'
43_LatOrig = -90
44_LatTile = _AtoZnoIO_.tillM
45_LonOrig = -180
46_LonTile = _AtoZnoIO_
47_M_      =  60000000000  # == 60_000_000_000 == 60 * pow(10, 9)
48_MaxPrec =  11
49_Radius_ =  Radius.__name__
50_Tile    =  15  # tile size in degrees
51
52_MaxLen  = _BaseLen + 2 * _MaxPrec
53_MinLen  = _BaseLen - 2
54
55_LatOrig_M_ = _LatOrig * _M_
56_LonOrig_M_ = _LonOrig * _M_
57
58_float_Tile   = _float(_Tile)
59_LatOrig_Tile = _float(_LatOrig) / _Tile
60_LonOrig_Tile = _float(_LonOrig) / _Tile
61
62
63def _2divmod3(x, Orig_M_):
64    '''(INTERNAL) Convert B{C{x}} to 3_tuple C{(tile, modulo, fraction)}/
65    '''
66    i = int(floor(x * _M_))
67    i, x = divmod(i - Orig_M_, _M_)
68    xt, xd = divmod(i, _Tile)
69    return xt, xd, x
70
71
72def _2fllh(lat, lon, height=None):
73    '''(INTERNAL) Convert lat, lon, height.
74    '''
75    # lat, lon = parseDMS2(lat, lon)
76    return (Lat(lat, Error=WGRSError),
77            Lon(lon, Error=WGRSError), height)
78
79
80def _2geostr2(georef):
81    '''(INTERNAL) Check a georef string.
82    '''
83    try:
84        n, geostr = len(georef), georef.upper()
85        p, o = divmod(n, 2)
86        if o or n < _MinLen or n > _MaxLen \
87             or geostr[:3] == _INV_ \
88             or not geostr.isalnum():
89            raise ValueError
90        return geostr, _2Precision(p - 1)
91
92    except (AttributeError, TypeError, ValueError) as x:
93        raise WGRSError(Georef.__name__, georef, txt=str(x))
94
95
96def _2Precision(precision):
97    '''(INTERNAL) Return a L{Precision_} instance.
98    '''
99    return Precision_(precision, Error=WGRSError, low=0, high=_MaxPrec)
100
101
102class WGRSError(_ValueError):
103    '''World Geographic Reference System (WGRS) encode, decode or other L{Georef} issue.
104    '''
105    pass
106
107
108class Georef(Str):
109    '''Georef class, a named C{str}.
110    '''
111    # no str.__init__ in Python 3
112    def __new__(cls, cll, precision=3, name=NN):
113        '''New L{Georef} from an other L{Georef} instance or georef
114           C{str} or from a C{LatLon} instance or lat-/longitude C{str}.
115
116           @arg cll: Cell or location (L{Georef} or C{str}, C{LatLon}
117                     or C{str}).
118           @kwarg precision: Optional, the desired georef resolution
119                             and length (C{int} 0..11), see function
120                             L{wgrs.encode} for more details.
121           @kwarg name: Optional name (C{str}).
122
123           @return: New L{Georef}.
124
125           @raise RangeError: Invalid B{C{cll}} lat- or longitude.
126
127           @raise TypeError: Invalid B{C{cll}}.
128
129           @raise WGRSError: INValid or non-alphanumeric B{C{cll}}.
130        '''
131        ll = p = None
132
133        if isinstance(cll, Georef):
134            g, p = _2geostr2(str(cll))
135
136        elif isstr(cll):
137            if ',' in cll:
138                lat, lon, h = _2fllh(*parse3llh(cll, height=None))
139                g  = encode(lat, lon, precision=precision, height=h)  # PYCHOK false
140                ll = lat, lon  # original lat, lon
141            else:
142                g = cll.upper()
143
144        else:  # assume LatLon
145            try:
146                lat, lon, h = _2fllh(cll.lat, cll.lon)
147                h = getattr(cll, _height_, h)
148            except AttributeError:
149                raise _xStrError(Georef, cll=cll)  # Error=WGRSError
150            g  = encode(lat, lon, precision=precision, height=h)  # PYCHOK false
151            ll = lat, lon  # original lat, lon
152
153        self = Str.__new__(cls, g, name=name or nameof(cll))
154        self._latlon    = ll
155        self._precision = p
156        return self
157
158    @Property_RO
159    def decoded3(self):
160        '''Get this georef's attributes (L{LatLonPrec3Tuple}).
161        '''
162        lat, lon = self.latlon
163        return LatLonPrec3Tuple(lat, lon, self.precision, name=self.name)
164
165    @Property_RO
166    def decoded5(self):
167        '''Get this georef's attributes (L{LatLonPrec5Tuple}) with
168           height and radius set to C{None} if missing.
169        '''
170        return self.decoded3.to5Tuple(self.height, self.radius)
171
172    @Property_RO
173    def _decoded5(self):
174        '''(INTERNAL) Initial L{LatLonPrec5Tuple}.
175        '''
176        return decode5(self)
177
178    @Property_RO
179    def height(self):
180        '''Get this georef's height in C{meter} or C{None} if missing.
181        '''
182        return self._decoded5.height
183
184    @Property_RO
185    def latlon(self):
186        '''Get this georef's (center) lat- and longitude (L{LatLon2Tuple}).
187        '''
188        lat, lon = self._latlon or self._decoded5[:2]
189        return LatLon2Tuple(lat, lon, name=self.name)
190
191    @Property_RO
192    def latlonheight(self):
193        '''Get this georef's (center) lat-, longitude and height (L{LatLon3Tuple}),
194           with height set to C{0} if missing.
195        '''
196        return self.latlon.to3Tuple(self.height or 0)
197
198    @Property_RO
199    def precision(self):
200        '''Get this georef's precision (C{int}).
201        '''
202        p = self._precision
203        return self._decoded5.precision if p is None else p
204
205    @Property_RO
206    def radius(self):
207        '''Get this georef's radius in C{meter} or C{None} if missing.
208        '''
209        return self._decoded5.radius
210
211    def toLatLon(self, LatLon=None, height=None, **LatLon_kwds):
212        '''Return (the center of) this georef cell as an instance
213           of the supplied C{LatLon} class.
214
215           @kwarg LatLon: Class to use (C{LatLon}) or C{None}.
216           @kwarg height: Optional height ({meter}).
217           @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
218                               arguments, ignored if C{B{LatLon} is None}.
219
220           @return: This georef location (B{C{LatLon}}) or if B{C{LatLon}}
221                    is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)}.
222
223           @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}}.
224        '''
225        if LatLon is None:
226            r = self.latlonheight if height is None else \
227                self.latlon.to3Tuple(height)
228        else:
229            h = (self.height or 0) if height is None else height
230            r =  LatLon(*self.latlon, height=h, **_xkwds(LatLon_kwds, name=self.name))
231        return r
232
233
234def decode3(georef, center=True):
235    '''Decode a C{georef} to lat-, longitude and precision.
236
237       @arg georef: To be decoded (L{Georef} or C{str}).
238       @kwarg center: If C{True} the center, otherwise the south-west,
239                      lower-left corner (C{bool}).
240
241       @return: A L{LatLonPrec3Tuple}C{(lat, lon, precision)}.
242
243       @raise WGRSError: Invalid B{C{georef}}, INValid, non-alphanumeric
244                         or odd length B{C{georef}}.
245    '''
246    def _digit(ll, g, i, m):
247        d = _Digits.find(g[i])
248        if d < 0 or d >= m:
249            raise _Error(i)
250        return ll * m + d
251
252    def _Error(i):
253        return WGRSError(Fmt.SQUARE(georef=i), georef)
254
255    def _index(chars, g, i):
256        k = chars.find(g[i])
257        if k < 0:
258            raise _Error(i)
259        return k
260
261    g, precision = _2geostr2(georef)
262    lon = _index(_LonTile, g, 0) + _LonOrig_Tile
263    lat = _index(_LatTile, g, 1) + _LatOrig_Tile
264
265    u = _1_0
266    if precision > 0:
267        lon = lon * _Tile + _index(_DegChar, g, 2)
268        lat = lat * _Tile + _index(_DegChar, g, 3)
269        m, p = 6, precision - 1
270        for i in range(_BaseLen, _BaseLen + p):
271            lon = _digit(lon, g, i,     m)
272            lat = _digit(lat, g, i + p, m)
273            u *= m
274            m  = _Base
275        u *= _Tile
276
277    if center:
278        lon = lon * _2_0 + _1_0
279        lat = lat * _2_0 + _1_0
280        u *= _2_0
281    u = _Tile / u
282    return LatLonPrec3Tuple(Lat(lat * u, Error=WGRSError),
283                            Lon(lon * u, Error=WGRSError),
284                            precision, name=nameof(georef))
285
286
287def decode5(georef, center=True):
288    '''Decode a C{georef} to lat-, longitude, precision, height and radius.
289
290       @arg georef: To be decoded (L{Georef} or C{str}).
291       @kwarg center: If C{True} the center, otherwise the south-west,
292                      lower-left corner (C{bool}).
293
294       @return: A L{LatLonPrec5Tuple}C{(lat, lon,
295                precision, height, radius)} where C{height} and/or
296                C{radius} are C{None} if missing.
297
298       @raise WGRSError: Invalid B{C{georef}}, INValid, non-alphanumeric
299                         or odd length B{C{georef}}.
300    '''
301    def _h2m(kft, name):
302        return Height(ft2m(kft * _1000_0), name=name, Error=WGRSError)
303
304    def _r2m(NM, name):
305        return Radius(NM / m2NM(1), name=name, Error=WGRSError)
306
307    def _split2(g, name, _2m):
308        i = max(g.find(name[0]), g.rfind(name[0]))
309        if i > _BaseLen:
310            return g[:i], _2m(int(g[i+1:]), _SPACE_(georef, name))
311        else:
312            return g, None
313
314    g = Str(georef, Error=WGRSError)
315
316    g, h = _split2(g, _Height_, _h2m)  # H is last
317    g, r = _split2(g, _Radius_, _r2m)  # R before H
318
319    return decode3(g, center=center).to5Tuple(h, r)
320
321
322def encode(lat, lon, precision=3, height=None, radius=None):  # MCCABE 14
323    '''Encode a lat-/longitude as a C{georef} of the given precision.
324
325       @arg lat: Latitude (C{degrees}).
326       @arg lon: Longitude (C{degrees}).
327       @kwarg precision: Optional, the desired C{georef} resolution and length
328                         (C{int} 0..11).
329       @kwarg height: Optional, height in C{meter}, see U{Designation of area
330                      <https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.
331       @kwarg radius: Optional, radius in C{meter}, see U{Designation of area
332                      <https://WikiPedia.org/wiki/World_Geographic_Reference_System>}.
333
334       @return: The C{georef} (C{str}).
335
336       @raise RangeError: Invalid B{C{lat}} or B{C{lon}}.
337
338       @raise WGRSError: Invalid B{C{precision}}, B{C{height}} or B{C{radius}}.
339
340       @note: The B{C{precision}} value differs from U{Georef<https://
341              GeographicLib.SourceForge.io/html/classGeographicLib_1_1Georef.html>}.
342              The C{georef} length is M{2 * (precision + 1)} and the
343              C{georef} resolution is I{15°} for B{C{precision}} 0, I{1°}
344              for 1, I{1′} for 2, I{0.1′} for 3, I{0.01′} for 4, ...
345              M{10**(2 - precision)}.
346    '''
347    def _option(name, m, m2_, K):
348        f = Scalar_(m, name=name, Error=WGRSError)
349        return NN(name[0].upper(), int(m2_(f * K) + _0_5))
350
351    p = _2Precision(precision)
352
353    lat, lon, _ = _2fllh(lat, lon)
354    if lat == _90_0:
355        lat *= EPS1_2
356
357    xt, xd, x = _2divmod3(lon, _LonOrig_M_)
358    yt, yd, y = _2divmod3(lat, _LatOrig_M_)
359
360    g = _LonTile[xt], _LatTile[yt]
361    if p > 0:
362        g += _DegChar[xd], _DegChar[yd]
363        p -= 1
364        if p > 0:
365            d =  pow(_Base, _MaxPrec - p)
366            x = _0wd(p, x // d)
367            y = _0wd(p, y // d)
368            g += x, y
369
370    if radius is not None:  # R before H
371        g += _option(_radius_, radius, m2NM, _1_0),
372    if height is not None:  # H is last
373        g += _option(_height_, height, m2ft, _0_001),
374
375    return NN.join(g)  # XXX Georef(''.join(g))
376
377
378def precision(res):
379    '''Determine the L{Georef} precision to meet a required (geographic)
380       resolution.
381
382       @arg res: The required resolution (C{degrees}).
383
384       @return: The L{Georef} precision (C{int} 0..11).
385
386       @raise ValueError: Invalid B{C{res}}.
387
388       @see: Function L{wgrs.encode} for more C{precision} details.
389    '''
390    r = Scalar_(res=res)
391    for p in range(_MaxPrec):
392        if resolution(p) <= r:
393            return p
394    return _MaxPrec
395
396
397def resolution(prec):
398    '''Determine the (geographic) resolution of a given L{Georef} precision.
399
400       @arg prec: The given precision (C{int}).
401
402       @return: The (geographic) resolution (C{degrees}).
403
404       @raise ValueError: Invalid B{C{prec}}.
405
406       @see: Function L{wgrs.encode} for more C{precision} details.
407    '''
408    p = Int(prec=prec, Error=WGRSError)
409    if p > 1:
410        r = _1_0 / (_60_0 * pow(_Base, min(p, _MaxPrec) - 1))
411    elif p < 1:
412        r = _float_Tile
413    else:
414        r = _1_0
415    return r
416
417
418__all__ += _ALL_OTHER(decode3, decode5,  # functions
419                      encode, precision, resolution)
420
421# **) MIT License
422#
423# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved.
424#
425# Permission is hereby granted, free of charge, to any person obtaining a
426# copy of this software and associated documentation files (the "Software"),
427# to deal in the Software without restriction, including without limitation
428# the rights to use, copy, modify, merge, publish, distribute, sublicense,
429# and/or sell copies of the Software, and to permit persons to whom the
430# Software is furnished to do so, subject to the following conditions:
431#
432# The above copyright notice and this permission notice shall be included
433# in all copies or substantial portions of the Software.
434#
435# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
436# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
437# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
438# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
439# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
440# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
441# OTHER DEALINGS IN THE SOFTWARE.
442