grib2io.utils

  1"""
  2Collection of utility functions to assist in the encoding and decoding
  3of GRIB2 Messages.
  4"""
  5
  6import datetime
  7import struct
  8from decimal import Decimal, localcontext
  9from typing import Dict, List, Optional, Tuple, Type, Union
 10
 11import numpy as np
 12from numpy.typing import ArrayLike
 13
 14from .. import tables
 15from .. import templates
 16
 17
 18def decimal_to_scaled_int(
 19    value: Union[float, str, int],
 20    scale_factor: Optional[int] = None,
 21) -> Tuple[int, int]:
 22    """
 23    Convert a float-like value to a scaled integer using the minimal decimal scaling factor.
 24
 25    The input value is internally converted to a `Decimal` to ensure precise scaling.
 26
 27    Parameters
 28    ----------
 29    value : float, str, or int
 30        The numeric value to scale.
 31    scaled_value : int
 32        The integer result of scaling the original value by `10**scale_factor`.
 33
 34    Returns
 35    -------
 36    scale_factor : int
 37        The smallest power of 10 such that `value * 10**scale_factor` is an exact integer.
 38    scaled_value : int
 39        The integer result of scaling the original value by `10**scale_factor`.
 40    """
 41    dec_value = Decimal(str(value))  # Preserve exact decimal representation
 42
 43    with localcontext() as ctx:
 44        ctx.prec = 28
 45
 46        if scale_factor is not None:
 47            scaled = dec_value * (10 ** scale_factor)
 48            if scaled != scaled.to_integral_value():
 49                raise ValueError(
 50                    f"Value {value} cannot be exactly scaled by 10^{scale_factor}"
 51                )
 52            return scale_factor, int(scaled)
 53        else:
 54            scale_factor = 0
 55            while dec_value != dec_value.to_integral_value():
 56                dec_value *= 10
 57                scale_factor += 1
 58                if scale_factor > 20:
 59                    raise ValueError(
 60                        f"Could not find exact scale factor for value {value} within bounds."
 61                    )
 62            return scale_factor, int(dec_value)
 63
 64
 65def int2bin(i: int, nbits: int=8, output: Union[Type[str], Type[List]]=str):
 66    """
 67    Convert integer to binary string or list
 68
 69    The struct module unpack using ">i" will unpack a 32-bit integer from a
 70    binary string.
 71
 72    Parameters
 73    ----------
 74    i
 75        Integer value to convert to binary representation.
 76    nbits : default=8
 77        Number of bits to return.  Valid values are 8 [DEFAULT], 16, 32, and
 78        64.
 79    output : default=str
 80        Return data as `str` [DEFAULT] or `list` (list of ints).
 81
 82    Returns
 83    -------
 84    int2bin
 85        `str` or `list` (list of ints) of binary representation of the integer
 86        value.
 87    """
 88    i = int(i) if not isinstance(i,int) else i
 89    assert nbits in [8,16,32,64]
 90    bitstr = "{0:b}".format(i).zfill(nbits)
 91    if output is str:
 92        return bitstr
 93    elif output is list:
 94        return [int(b) for b in bitstr]
 95
 96
 97def ieee_float_to_int(f):
 98    """
 99    Convert an IEEE 754 32-bit float to a 32-bit integer.
100
101    Parameters
102    ----------
103    f : float
104        Floating-point value.
105
106    Returns
107    -------
108    ieee_float_to_int
109        `numpy.int32` representation of an IEEE 32-bit float.
110    """
111    i = struct.unpack('>i',struct.pack('>f',np.float32(f)))[0]
112    return np.int32(i)
113
114
115def ieee_int_to_float(i):
116    """
117    Convert a 32-bit integer to an IEEE 32-bit float.
118
119    Parameters
120    ----------
121    i : int
122        Integer value.
123
124    Returns
125    -------
126    ieee_int_to_float
127        `numpy.float32` representation of a 32-bit int.
128    """
129    f = struct.unpack('>f',struct.pack('>i',np.int32(i)))[0]
130    return np.float32(f)
131
132
133def get_leadtime(pdtn: int, pdt: ArrayLike) -> datetime.timedelta:
134    """
135    Compute lead time as a datetime.timedelta object.
136
137    Using information from GRIB2 Product Definition Template
138    Number, and Product Definition Template (Section 4).
139
140    Parameters
141    ----------
142    pdtn
143        GRIB2 Product Definition Template Number
144    pdt
145        Sequence containing GRIB2 Product Definition Template (Section 4).
146
147    Returns
148    -------
149    leadTime
150        datetime.timedelta object representing the lead time of the GRIB2 message.
151    """
152    lt = tables.get_value_from_table(pdt[templates.UnitOfForecastTime._key[pdtn]], 'scale_time_seconds')
153    lt *= pdt[templates.ValueOfForecastTime._key[pdtn]]
154    return datetime.timedelta(seconds=int(lt))
155
156
157def get_duration(pdtn: int, pdt: ArrayLike) -> datetime.timedelta:
158    """
159    Compute a time duration as a datetime.timedelta.
160
161    Uses information from Product Definition Template Number, and Product
162    Definition Template (Section 4).
163
164    Parameters
165    ----------
166    pdtn
167        GRIB2 Product Definition Template Number
168    pdt
169        Sequence containing GRIB2 Product Definition Template (Section 4).
170
171    Returns
172    -------
173    get_duration
174        datetime.timedelta object representing the time duration of the GRIB2
175        message.
176    """
177    if pdtn in templates._timeinterval_pdtns:
178        ntime = pdt[templates.NumberOfTimeRanges._key[pdtn]]
179        duration_unit = tables.get_value_from_table(
180            pdt[templates.UnitOfTimeRangeOfStatisticalProcess._key[pdtn]],
181            'scale_time_seconds')
182        d = ntime * duration_unit * pdt[
183            templates.TimeRangeOfStatisticalProcess._key[pdtn]]
184    else:
185        d = 0
186    return datetime.timedelta(seconds=int(d))
187
188
189def decode_wx_strings(lus: bytes) -> Dict[int, str]:
190    """
191    Decode GRIB2 Local Use Section to obtain NDFD/MDL Weather Strings.
192
193    The decode procedure is defined
194    [here](https://vlab.noaa.gov/web/mdl/nbm-gmos-grib2-wx-info).
195
196    Parameters
197    ----------
198    lus
199        GRIB2 Local Use Section containing NDFD weather strings.
200
201    Returns
202    -------
203    decode_wx_strings
204        Dict of NDFD/MDL weather strings. Keys are an integer value that
205        represent the sequential order of the key in the packed local use
206        section and the value is the weather key.
207    """
208    assert lus[0] == 1
209    # Unpack information related to the simple packing method
210    # the packed weather string data.
211    ngroups = struct.unpack('>H',lus[1:3])[0]
212    nvalues = struct.unpack('>i',lus[3:7])[0]
213    refvalue = struct.unpack('>i',lus[7:11])[0]
214    dsf = struct.unpack('>h',lus[11:13])[0]
215    nbits = lus[13]
216    datatype = lus[14]
217    if datatype == 0: # Floating point
218        refvalue = np.float32(ieee_int_to_float(refvalue)*10**-dsf)
219    elif datatype == 1: # Integer
220        refvalue = np.int32(ieee_int_to_float(refvalue)*10**-dsf)
221    # Upack each byte starting at byte 15 to end of the local use
222    # section, create a binary string and append to the full
223    # binary string.
224    b = ''
225    for i in range(15,len(lus)):
226        iword = struct.unpack('>B',lus[i:i+1])[0]
227        b += bin(iword).split('b')[1].zfill(8)
228    # Iterate over the binary string (b). For each nbits
229    # chunk, convert to an integer, including the refvalue,
230    # and then convert the int to an ASCII character, then
231    # concatenate to wxstring.
232    wxstring = ''
233    for i in range(0,len(b),nbits):
234        wxstring += chr(int(b[i:i+nbits],2)+refvalue)
235    # Return string as list, split by null character.
236    #return list(filter(None,wxstring.split('\0')))
237    return {n:k for n,k in enumerate(list(filter(None,wxstring.split('\0'))))}
238
239
240def get_wgrib2_prob_string(
241    probtype: int,
242    sfacl: int,
243    svall: int,
244    sfacu: int,
245    svalu: int,
246) -> str:
247    """
248    Return a wgrib2-styled string of probabilistic threshold information.
249
250    Logic from wgrib2 source,
251    [Prob.c](https://github.com/NOAA-EMC/NCEPLIBS-wgrib2/blob/develop/wgrib2/Prob.c),
252    is replicated here.
253
254    Parameters
255    ----------
256    probtype
257        Type of probability (Code Table 4.9).
258    sfacl
259        Scale factor of lower limit.
260    svall
261        Scaled value of lower limit.
262    sfacu
263        Scale factor of upper limit.
264    svalu
265        Scaled value of upper limit.
266
267    Returns
268    -------
269    get_wgrib2_prob_string
270        wgrib2-formatted string of probability threshold.
271    """
272    probstr = ''
273    if sfacl == -127: sfacl = 0
274    if sfacu == -127: sfacu = 0
275    lower = svall/(10**sfacl)
276    upper = svalu/(10**sfacu)
277    if probtype == 0:
278        probstr = 'prob <%g' % (lower)
279    elif probtype == 1:
280        probstr = 'prob >%g' % (upper)
281    elif probtype == 2:
282        if lower == upper:
283            probstr = 'prob =%g' % (lower)
284        else:
285            probstr = 'prob >=%g <%g' % (lower,upper)
286    elif probtype == 3:
287        probstr = 'prob >%g' % (lower)
288    elif probtype == 4:
289        probstr = 'prob <%g' % (upper)
290    else:
291        probstr = ''
292    return probstr