86class open():
 87    """
 88    GRIB2 File Object.
 89
 90    A physical file can contain one or more GRIB2 messages.  When instantiated,
 91    class `grib2io.open`, the file named `filename` is opened for reading (`mode
 92    = 'r'`) and is automatically indexed.  The indexing procedure reads some of
 93    the GRIB2 metadata for all GRIB2 Messages.  A GRIB2 Message may contain
 94    submessages whereby Section 2-7 can be repeated.  grib2io accommodates for
 95    this by flattening any GRIB2 submessages into multiple individual messages.
 96
 97    It is important to note that GRIB2 files from some Meteorological agencies
 98    contain other data than GRIB2 messages.  GRIB2 files from ECMWF can contain
 99    GRIB1 and GRIB2 messages.  grib2io checks for these and safely ignores them.
100
101    Attributes
102    ----------
103    closed : bool
104        `True` is file handle is close; `False` otherwise.
105    current_message : int
106        Current position of the file in units of GRIB2 Messages.
107    levels : tuple
108        Tuple containing a unique list of wgrib2-formatted level/layer strings.
109    messages : int
110        Count of GRIB2 Messages contained in the file.
111    mode : str
112        File IO mode of opening the file.
113    name : str
114        Full path name of the GRIB2 file.
115    size : int
116        Size of the file in units of bytes.
117    variables : tuple
118        Tuple containing a unique list of variable short names (i.e. GRIB2
119        abbreviation names).
120    """
121
122    __slots__ = ('_fileid', '_filehandle', '_hasindex', '_index', '_nodata',
123                 '_pos', 'closed', 'current_message', 'messages', 'mode',
124                 'name', 'size')
125
126    def __init__(
127        self,
128        filename: Union[bytes, str, Path],
129        mode: Literal["r", "w", "x"] = "r",
130        **kwargs,
131        ):
132        """
133        Initialize GRIB2 File object instance.
134
135        Parameters
136        ----------
137        filename
138            File name containing GRIB2 messages.
139            OR
140            bytes
141        mode: default="r"
142            File access mode where "r" opens the files for reading only; "w"
143            opens the file for overwriting and "x" for writing to a new file.
144        """
145
146        # Manage keywords
147        if "_xarray_backend" not in kwargs:
148            kwargs["_xarray_backend"] = False
149            self._nodata = False
150        else:
151            self._nodata = kwargs["_xarray_backend"]
152
153        # All write modes are read/write.
154        # All modes are binary.
155        if mode in ("a", "x", "w"):
156            mode += "+"
157        mode = mode + "b"
158
159        if isinstance(filename, bytes):
160            if filename[:2] == _GZIP_HEADER:
161                filename = gzip.decompress(filename)
162            if filename[:4]+filename[-4:] != b'GRIB7777':
163                raise ValueError("Invalid GRIB bytes")
164            self._filehandle = BytesIO(filename)
165            self.name = "<in-memory-file>"
166            self.size = len(filename)
167            self._fileid = hashlib.sha1((self.name+str(self.size)).encode('ASCII')).hexdigest()
168
169        # Some GRIB2 files are gzipped, so check for that here, but
170        # raise error when using xarray backend.
171        else:
172            if 'r' in mode:
173                self._filehandle = builtins.open(filename, mode=mode)
174                # Gzip files contain a 2-byte header b'\x1f\x8b'.
175                if self._filehandle.read(2) == _GZIP_HEADER:
176                    self._filehandle.close()
177                    if kwargs["_xarray_backend"]:
178                        raise RuntimeError('Gzip GRIB2 files are not supported by the Xarray backend.')
179                    self._filehandle = gzip.open(filename, mode=mode)
180                else:
181                    self._filehandle = builtins.open(filename, mode=mode)
182            else:
183                self._filehandle = builtins.open(filename, mode=mode)
184                self.name = os.path.abspath(filename)
185            self.name = os.path.abspath(filename)
186            fstat = os.stat(self.name)
187            self.size = fstat.st_size
188            self._fileid = hashlib.sha1((self.name+str(fstat.st_ino)+
189                                     str(self.size)).encode('ASCII')).hexdigest()
190
191        self._hasindex = False
192        self._index = {}
193        self.mode = mode
194        self.messages = 0
195        self.current_message = 0
196        self.closed = self._filehandle.closed
197        if 'r' in self.mode:
198            self._build_index()
199        # FIX: Cannot perform reads on mode='a'
200        #if 'a' in self.mode and self.size > 0: self._build_index()
201
202
203    def __delete__(self, instance):
204        self.close()
205        del self._index
206
207
208    def __enter__(self):
209        return self
210
211
212    def __exit__(self, atype, value, traceback):
213        self.close()
214
215
216    def __iter__(self):
217        yield from self._index['msg']
218
219
220    def __len__(self):
221        return self.messages
222
223
224    def __repr__(self):
225        strings = []
226        for k in self.__slots__:
227            if k.startswith('_'): continue
228            strings.append('%s = %s\n'%(k,eval('self.'+k)))
229        return ''.join(strings)
230
231
232    def __getitem__(self, key):
233        if isinstance(key,int):
234            if abs(key) >= len(self._index['msg']):
235                raise IndexError("index out of range")
236            else:
237                return self._index['msg'][key]
238        elif isinstance(key,str):
239            return self.select(shortName=key)
240        elif isinstance(key,slice):
241            return self._index['msg'][key]
242        else:
243            raise KeyError('Key must be an integer, slice, or GRIB2 variable shortName.')
244
245
246    def _build_index(self):
247        """Perform indexing of GRIB2 Messages."""
248        # Initialize index dictionary
249        if not self._hasindex:
250            self._index['sectionOffset'] = []
251            self._index['sectionSize'] = []
252            self._index['msgSize'] = []
253            self._index['msgNumber'] = []
254            self._index['msg'] = []
255            self._index['isSubmessage'] = []
256            self._hasindex = True
257
258        # Iterate
259        while True:
260            try:
261                # Initialize
262                bmapflag = None
263                pos = self._filehandle.tell()
264                section2 = b''
265                trailer = b''
266                _secpos = dict.fromkeys(range(8))
267                _secsize = dict.fromkeys(range(8))
268                _isSubmessage = False
269
270                # Ignore headers (usually text) that are not part of the GRIB2
271                # file.  For example, NAVGEM files have a http header at the
272                # beginning that needs to be ignored.
273
274                # Read a byte at a time until "GRIB" is found.  Using
275                # "wgrib2" on a NAVGEM file, the header was 421 bytes and
276                # decided to go to 2048 bytes to be safe. For normal GRIB2
277                # files this should be quick and break out of the first
278                # loop when "test_offset" is 0.
279                for test_offset in range(2048):
280                    self._filehandle.seek(pos + test_offset)
281                    header = struct.unpack(">I", self._filehandle.read(4))[0]
282                    if header.to_bytes(4, "big") == b"GRIB":
283                        pos = pos + test_offset
284                        break
285                else:
286                    # NOTE: Coming here means that no "GRIB" message identifier
287                    # was found in the previous 2048 bytes. So here we continue
288                    # the while True loop.
289                    continue
290
291                # Read the rest of Section 0 using struct.
292                _secpos[0] = self._filehandle.tell()-4
293                _secsize[0] = 16
294                secmsg = self._filehandle.read(12)
295                section0 = np.concatenate(([header],list(struct.unpack('>HBBQ',secmsg))),dtype=np.int64)
296
297                # Make sure message is GRIB2.
298                if section0[3] != 2:
299                    # Check for GRIB1 and ignore.
300                    if secmsg[3] == 1:
301                        warnings.warn("GRIB version 1 message detected.  Ignoring...")
302                        grib1_size = int.from_bytes(secmsg[:3],"big")
303                        self._filehandle.seek(self._filehandle.tell()+grib1_size-16)
304                        continue
305                    else:
306                        raise ValueError("Bad GRIB version number.")
307
308                # Read and unpack sections 1 through 8 which all follow a
309                # pattern of section size, number, and content.
310                while 1:
311                    # Read first 5 bytes of the section which contains the size
312                    # of the section (4 bytes) and the section number (1 byte).
313                    secmsg = self._filehandle.read(5)
314                    secsize, secnum = struct.unpack('>iB',secmsg)
315
316                    # Record the offset of the section number and "append" the
317                    # rest of the section to secmsg.
318                    _secpos[secnum] = self._filehandle.tell()-5
319                    _secsize[secnum] = secsize
320                    if secnum in {1,3,4,5}:
321                        secmsg += self._filehandle.read(secsize-5)
322                    grbpos = 0
323
324                    # Unpack section
325                    if secnum == 1:
326                        # Unpack Section 1
327                        section1, grbpos = g2clib.unpack1(secmsg)
328                    elif secnum == 2:
329                        # Unpack Section 2
330                        section2 = self._filehandle.read(secsize-5)
331                    elif secnum == 3:
332                        # Unpack Section 3
333                        gds, gdt, deflist, grbpos = g2clib.unpack3(secmsg)
334                        gds = gds.tolist()
335                        gdt = gdt.tolist()
336                        section3 = np.concatenate((gds,gdt))
337                        section3 = np.where(section3==4294967295,-1,section3)
338                    elif secnum == 4:
339                        # Unpack Section 4
340                        pdtnum, pdt, coordlist, numcoord, grbpos = g2clib.unpack4(secmsg)
341                        pdt = pdt.tolist()
342                        section4 = np.concatenate((np.array((numcoord,pdtnum)),pdt))
343                    elif secnum == 5:
344                        # Unpack Section 5
345                        drtn, drt, npts, self._pos = g2clib.unpack5(secmsg)
346                        section5 = np.concatenate((np.array((npts,drtn)),drt))
347                        section5 = np.where(section5==4294967295,-1,section5)
348                    elif secnum == 6:
349                        # Unpack Section 6 - Just the bitmap flag
350                        bmapflag = struct.unpack('>B',self._filehandle.read(1))[0]
351                        if bmapflag == 0:
352                            bmappos = self._filehandle.tell()-6
353                        elif bmapflag == 254:
354                            # Do this to keep the previous position value
355                            pass
356                        else:
357                            bmappos = None
358                        self._filehandle.seek(self._filehandle.tell()+secsize-6)
359                    elif secnum == 7:
360                        # Do not unpack section 7 here, but move the file pointer
361                        # to after section 7.
362                        self._filehandle.seek(self._filehandle.tell()+secsize-5)
363
364                        # Update the file index.
365                        self.messages += 1
366                        self._index['sectionOffset'].append(copy.deepcopy(_secpos))
367                        self._index['sectionSize'].append(copy.deepcopy(_secsize))
368                        self._index['msgSize'].append(section0[-1])
369                        self._index['msgNumber'].append(self.messages)
370                        self._index['isSubmessage'].append(_isSubmessage)
371
372                        # Create Grib2Message with data.
373                        msg = Grib2Message(section0,section1,section2,section3,section4,section5,bmapflag)
374                        msg._msgnum = self.messages-1
375                        msg._deflist = deflist
376                        msg._coordlist = coordlist
377                        if not self._nodata:
378                            msg._ondiskarray = Grib2MessageOnDiskArray((msg.ny,msg.nx), 2,
379                                                                TYPE_OF_VALUES_DTYPE[msg.typeOfValues],
380                                                                self._filehandle,
381                                                                msg, pos, _secpos[6], _secpos[7])
382                        self._index['msg'].append(msg)
383
384                        # If here, then we have moved through GRIB2 section 1-7.
385                        # Now we need to check the next 4 bytes after section 7.
386                        trailer = struct.unpack('>i',self._filehandle.read(4))[0]
387
388                        # If we reach the GRIB2 trailer string ('7777'), then we
389                        # can break begin processing the next GRIB2 message.  If
390                        # not, then we continue within the same iteration to
391                        # process a GRIB2 submessage.
392                        if trailer.to_bytes(4, "big") == b'7777':
393                            break
394                        else:
395                            # If here, trailer should be the size of the first
396                            # section of the next submessage, then the next byte
397                            # is the section number.  Check this value.
398                            nextsec = struct.unpack('>B',self._filehandle.read(1))[0]
399                            if nextsec not in {2,3,4}:
400                                raise ValueError("Bad GRIB2 message structure.")
401                            self._filehandle.seek(self._filehandle.tell()-5)
402                            _isSubmessage = True
403                            continue
404                    else:
405                        raise ValueError("Bad GRIB2 section number.")
406
407            except(struct.error):
408                if 'r' in self.mode:
409                    self._filehandle.seek(0)
410                break
411
412
413    @property
414    def levels(self):
415        """Provides a unique tuple of level strings."""
416        if self._hasindex and not self._nodata:
417            return tuple(sorted(set([msg.level for msg in self._index['msg']])))
418        else:
419            return None
420
421
422    @property
423    def variables(self):
424        """Provides a unique tuple of variable shortName strings."""
425        if self._hasindex and not self._nodata:
426            return tuple(sorted(set([msg.shortName for msg in self._index['msg']])))
427        else:
428            return None
429
430
431    def close(self):
432        """Close the file handle."""
433        if not self._filehandle.closed:
434            self.messages = 0
435            self.current_message = 0
436            self._filehandle.close()
437            self.closed = self._filehandle.closed
438
439
440    def read(self, size: Optional[int]=None):
441        """
442        Read size amount of GRIB2 messages from the current position.
443
444        If no argument is given, then size is None and all messages are returned
445        from the current position in the file. This read method follows the
446        behavior of Python's builtin open() function, but whereas that operates
447        on units of bytes, we operate on units of GRIB2 messages.
448
449        Parameters
450        ----------
451        size: default=None
452            The number of GRIB2 messages to read from the current position. If
453            no argument is give, the default value is None and remainder of
454            the file is read.
455
456        Returns
457        -------
458        read
459            ``Grib2Message`` object when size = 1 or a list of Grib2Messages
460            when size > 1.
461        """
462        if size is not None and size < 0:
463            size = None
464        if size is None or size > 1:
465            start = self.tell()
466            stop = self.messages if size is None else start+size
467            if size is None:
468                self.current_message = self.messages-1
469            else:
470                self.current_message += size
471            return self._index['msg'][slice(start,stop,1)]
472        elif size == 1:
473            self.current_message += 1
474            return self._index['msg'][self.current_message]
475        else:
476            None
477
478
479    def seek(self, pos: int):
480        """
481        Set the position within the file in units of GRIB2 messages.
482
483        Parameters
484        ----------
485        pos
486            The GRIB2 Message number to set the file pointer to.
487        """
488        if self._hasindex:
489            self._filehandle.seek(self._index['sectionOffset'][0][pos])
490            self.current_message = pos
491
492
493    def tell(self):
494        """Returns the position of the file in units of GRIB2 Messages."""
495        return self.current_message
496
497
498    def select(self, **kwargs):
499        """Select GRIB2 messages by `Grib2Message` attributes."""
500        # TODO: Added ability to process multiple values for each keyword (attribute)
501        idxs = []
502        nkeys = len(kwargs.keys())
503        for k,v in kwargs.items():
504            for m in self._index['msg']:
505                if hasattr(m,k) and getattr(m,k) == v: idxs.append(m._msgnum)
506        idxs = np.array(idxs,dtype=np.int32)
507        return [self._index['msg'][i] for i in [ii[0] for ii in collections.Counter(idxs).most_common() if ii[1] == nkeys]]
508
509
510    def write(self, msg):
511        """
512        Writes GRIB2 message object to file.
513
514        Parameters
515        ----------
516        msg
517            GRIB2 message objects to write to file.
518        """
519        if isinstance(msg,list):
520            for m in msg:
521                self.write(m)
522            return
523
524        if issubclass(msg.__class__,_Grib2Message):
525            if hasattr(msg,'_msg'):
526                self._filehandle.write(msg._msg)
527            else:
528                if msg._signature != msg._generate_signature():
529                    msg.pack()
530                    self._filehandle.write(msg._msg)
531                else:
532                    if hasattr(msg._data,'filehandle'):
533                        msg._data.filehandle.seek(msg._data.offset)
534                        self._filehandle.write(msg._data.filehandle.read(msg.section0[-1]))
535                    else:
536                        msg.pack()
537                        self._filehandle.write(msg._msg)
538            self.flush()
539            self.size = os.path.getsize(self.name)
540            self._filehandle.seek(self.size-msg.section0[-1])
541            self._build_index()
542        else:
543            raise TypeError("msg must be a Grib2Message object.")
544        return
545
546
547    def flush(self):
548        """Flush the file object buffer."""
549        self._filehandle.flush()
550
551
552    def levels_by_var(self, name: str):
553        """
554        Return a list of level strings given a variable shortName.
555
556        Parameters
557        ----------
558        name
559            Grib2Message variable shortName
560
561        Returns
562        -------
563        levels_by_var
564            A list of unique level strings.
565        """
566        return list(sorted(set([msg.level for msg in self.select(shortName=name)])))
567
568
569    def vars_by_level(self, level: str):
570        """
571        Return a list of variable shortName strings given a level.
572
573        Parameters
574        ----------
575        level
576            Grib2Message variable level
577
578        Returns
579        -------
580        vars_by_level
581            A list of unique variable shortName strings.
582        """
583        return list(sorted(set([msg.shortName for msg in self.select(level=level)])))

34def show_config():
35    """Print grib2io build configuration information."""
36    print(f'grib2io version {__version__} Configuration:')
37    print(f'')
38    print(f'NCEPLIBS-g2c library version: {__g2clib_version__}')
39    print(f'\tStatic library: {g2c_static}')
40    print(f'\tJPEG compression support: {has_jpeg_support}')
41    print(f'\tPNG compression support: {has_png_support}')
42    print(f'\tAEC compression support: {has_aec_support}')
43    print(f'')
44    print(f'NCEPLIPS-ip support: {has_interpolation}')
45    print(f'\tStatic library: {ip_static}')
46    print(f'\tOpenMP support: {has_openmp_support}')
47    print(f'')
48    print(f'Static libs:')
49    for lib in extra_objects:
50        print(f'\t{lib}')
51    print(f'')
52    print(f'NCEP GRIB2 Table Version: {_ncep_grib2_table_version}')

1698def interpolate(a,
1699    method: Union[int, str],
1700    grid_def_in,
1701    grid_def_out,
1702    method_options=None,
1703    num_threads=1,
1704):
1705    """
1706    This is the module-level interpolation function.
1707
1708    This interfaces with the 4-byte (32-bit float) [NCEPLIBS-ip library](https://github.com/NOAA-EMC/NCEPLIBS-ip)
1709    through grib2io's internal iplib Cython extension module.
1710
1711    Parameters
1712    ----------
1713    a : numpy.ndarray or tuple
1714        Input data.  If `a` is a `numpy.ndarray`, scalar interpolation will be
1715        performed.  If `a` is a `tuple`, then vector interpolation will be
1716        performed with the assumption that u = a[0] and v = a[1] and are both
1717        `numpy.ndarray`.
1718
1719        These data are expected to be in 2-dimensional form with shape (ny, nx)
1720        or 3-dimensional (:, ny, nx) where the 1st dimension represents another
1721        spatial, temporal, or classification (i.e. ensemble members) dimension.
1722        The function will properly flatten the (ny,nx) dimensions into (nx * ny)
1723        acceptable for input into the interpolation subroutines. If needed, these
1724        data will be converted to `np.float32`.
1725    method
1726        Interpolate method to use. This can either be an integer or string using
1727        the following mapping:
1728
1729        | Interpolate Scheme | Integer Value |
1730        | :---:              | :---:         |
1731        | 'bilinear'         | 0             |
1732        | 'bicubic'          | 1             |
1733        | 'neighbor'         | 2             |
1734        | 'budget'           | 3             |
1735        | 'spectral'         | 4             |
1736        | 'neighbor-budget'  | 6             |
1737
1738    grid_def_in : grib2io.Grib2GridDef
1739        Grib2GridDef object for the input grid.
1740    grid_def_out : grib2io.Grib2GridDef
1741        Grib2GridDef object for the output grid or station points.
1742    method_options : list of ints, optional
1743        Interpolation options. See the NCEPLIBS-ip documentation for
1744        more information on how these are used.
1745    num_threads : int, optional
1746        Number of OpenMP threads to use for interpolation. The default
1747        value is 1. If NCEPLIBS-ip and grib2io's iplib extension module
1748        was not built with OpenMP, then this keyword argument and value
1749        will have no impact.
1750
1751    Returns
1752    -------
1753    interpolate
1754        Returns a `numpy.ndarray` of dtype `np.float32` when scalar interpolation
1755        is performed or a `tuple` of `numpy.ndarray`s when vector interpolation is
1756        performed with the assumptions that 0-index is the interpolated u and
1757        1-index is the interpolated v.
1758    """
1759
1760    from . import iplib
1761
1762    prev_num_threads = 1
1763    try:
1764        prev_num_threads = iplib.openmp_get_num_threads()
1765        iplib.openmp_set_num_threads(num_threads)
1766    except(AttributeError):
1767        pass
1768
1769    print(f"grib2io.interpolate thread report: OpenMP num threads = {iplib.openmp_get_num_threads()}")
1770
1771    if isinstance(method,int) and method not in _interp_schemes.values():
1772        raise ValueError('Invalid interpolation method.')
1773    elif isinstance(method,str):
1774        if method in _interp_schemes.keys():
1775            method = _interp_schemes[method]
1776        else:
1777            raise ValueError('Invalid interpolation method.')
1778
1779    if method_options is None:
1780        method_options = np.zeros((20),dtype=np.int32)
1781        if method in {3,6}:
1782            method_options[0:2] = -1
1783
1784    mi = grid_def_in.npoints
1785    mo = grid_def_out.npoints
1786
1787    # Adjust shape of input array(s)
1788    a, newshp = _adjust_array_shape_for_interp(a, grid_def_in, grid_def_out)
1789
1790    # Call interpolation subroutines according to type of a.
1791    if isinstance(a,np.ndarray):
1792        # Scalar
1793        km = a.shape[0]
1794        if np.any(np.isnan(a)):
1795            ibi = np.ones((km), dtype=np.int32)
1796            li = np.where(np.isnan(a),0,1).astype(np.uint8)
1797        else:
1798            ibi = np.zeros((km), dtype=np.int32)
1799            li = np.zeros(a.shape,dtype=np.uint8)
1800        no, rlat, rlon, ibo, lo, go, iret = iplib.interpolate_scalar(
1801            method,
1802            method_options,
1803            grid_def_in.gdtn,
1804            np.array(grid_def_in.gdt, dtype=np.int32),
1805            grid_def_out.gdtn,
1806            np.array(grid_def_out.gdt, dtype=np.int32),
1807            mi,
1808            mo,
1809            km,
1810            ibi,
1811            li,
1812            a,
1813        )
1814        out = np.where(lo==0,np.nan,go).reshape(newshp)
1815    elif isinstance(a,tuple):
1816        # Vector
1817        km = a[0].shape[0]
1818        if np.any(np.isnan(a)):
1819            ibi = np.ones((km), dtype=np.int32)
1820            li = np.where(np.isnan(a),0,1).astype(np.uint8)
1821        else:
1822            ibi = np.zeros((km), dtype=np.int32)
1823            li = np.zeros(a[0].shape,dtype=np.uint8)
1824        no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret = iplib.interpolate_vector(
1825            method,
1826            method_options,
1827            grid_def_in.gdtn,
1828            np.array(grid_def_in.gdt, dtype=np.int32),
1829            grid_def_out.gdtn,
1830            np.array(grid_def_out.gdt, dtype=np.int32),
1831            mi,
1832            mo,
1833            km,
1834            ibi,
1835            li,
1836            a[0],
1837            a[1],
1838        )
1839        uo = np.where(lo==0,np.nan,uo).reshape(newshp)
1840        vo = np.where(lo==0,np.nan,vo).reshape(newshp)
1841        out = (uo,vo)
1842
1843    try:
1844        iplib.openmp_set_num_threads(prev_num_threads)
1845    except(AttributeError):
1846        pass
1847
1848    return out

1851def interpolate_to_stations(
1852    a,
1853    method: Union[int, str],
1854    grid_def_in,
1855    lats,
1856    lons,
1857    method_options=None,
1858    num_threads=1
1859):
1860    """
1861    Module-level interpolation function for interpolation to stations.
1862
1863    Interfaces with the 4-byte (32-bit float) [NCEPLIBS-ip library](https://github.com/NOAA-EMC/NCEPLIBS-ip)
1864    via grib2io's iplib Cython exntension module. It supports scalar and
1865    vector interpolation according to the type of object a.
1866
1867    Parameters
1868    ----------
1869    a : numpy.ndarray or tuple
1870        Input data.  If `a` is a `numpy.ndarray`, scalar interpolation will be
1871        performed.  If `a` is a `tuple`, then vector interpolation will be
1872        performed with the assumption that u = a[0] and v = a[1] and are both
1873        `numpy.ndarray`.
1874
1875        These data are expected to be in 2-dimensional form with shape (ny, nx)
1876        or 3-dimensional (:, ny, nx) where the 1st dimension represents another
1877        spatial, temporal, or classification (i.e. ensemble members) dimension.
1878        The function will properly flatten the (ny,nx) dimensions into (nx * ny)
1879        acceptable for input into the interpolation subroutines. If needed, these
1880        data will be converted to `np.float32`.
1881    method
1882        Interpolate method to use. This can either be an integer or string using
1883        the following mapping:
1884
1885        | Interpolate Scheme | Integer Value |
1886        | :---:              | :---:         |
1887        | 'bilinear'         | 0             |
1888        | 'bicubic'          | 1             |
1889        | 'neighbor'         | 2             |
1890        | 'budget'           | 3             |
1891        | 'spectral'         | 4             |
1892        | 'neighbor-budget'  | 6             |
1893
1894    grid_def_in : grib2io.Grib2GridDef
1895        Grib2GridDef object for the input grid.
1896    lats : numpy.ndarray or list
1897        Latitudes for station points
1898    lons : numpy.ndarray or list
1899        Longitudes for station points
1900    method_options : list of ints, optional
1901        Interpolation options. See the NCEPLIBS-ip documentation for
1902        more information on how these are used.
1903    num_threads : int, optional
1904        Number of OpenMP threads to use for interpolation. The default
1905        value is 1. If NCEPLIBS-ip and grib2io's iplib extension module
1906        was not built with OpenMP, then this keyword argument and value
1907        will have no impact.
1908
1909    Returns
1910    -------
1911    interpolate_to_stations
1912        Returns a `numpy.ndarray` of dtype `np.float32` when scalar
1913        interpolation is performed or a `tuple` of `numpy.ndarray`s
1914        when vector interpolation is performed with the assumptions
1915        that 0-index is the interpolated u and 1-index is the
1916        interpolated v.
1917    """
1918    from . import iplib
1919
1920    prev_num_threads = 1
1921    try:
1922        prev_num_threads = iplib.openmp_get_num_threads()
1923        iplib.openmp_set_num_threads(num_threads)
1924    except(AttributeError):
1925        pass
1926
1927    if isinstance(method,int) and method not in _interp_schemes.values():
1928        raise ValueError('Invalid interpolation method.')
1929    elif isinstance(method,str):
1930        if method in _interp_schemes.keys():
1931            method = _interp_schemes[method]
1932        else:
1933            raise ValueError('Invalid interpolation method.')
1934
1935    if method_options is None:
1936        method_options = np.zeros((20),dtype=np.int32)
1937        if method in {3,6}:
1938            method_options[0:2] = -1
1939
1940    # Check lats and lons
1941    if isinstance(lats,list):
1942        nlats = len(lats)
1943    elif isinstance(lats,np.ndarray) and len(lats.shape) == 1:
1944        nlats = lats.shape[0]
1945    else:
1946        raise ValueError('Station latitudes must be a list or 1-D NumPy array.')
1947    if isinstance(lons,list):
1948        nlons = len(lons)
1949    elif isinstance(lons,np.ndarray) and len(lons.shape) == 1:
1950        nlons = lons.shape[0]
1951    else:
1952        raise ValueError('Station longitudes must be a list or 1-D NumPy array.')
1953    if nlats != nlons:
1954        raise ValueError('Station lats and lons must be same size.')
1955
1956    mi = grid_def_in.npoints
1957    mo = nlats
1958
1959    # Adjust shape of input array(s)
1960    a, newshp = _adjust_array_shape_for_interp_stations(a, grid_def_in, mo)
1961
1962    # Use gdtn = -1 for stations and an empty template array
1963    gdtn = -1
1964    gdt = np.zeros((200),dtype=np.int32)
1965
1966    # Call interpolation subroutines according to type of a.
1967    if isinstance(a,np.ndarray):
1968        # Scalar
1969        km = a.shape[0]
1970        ibi = np.zeros((km), dtype=np.int32)
1971        li = np.zeros(a.shape,dtype=np.uint8)
1972        no, rlat, rlon, ibo, lo, go, iret = iplib.interpolate_scalar(
1973            method,
1974            method_options,
1975            grid_def_in.gdtn,
1976            np.array(grid_def_in.gdt, dtype=np.int32),
1977            gdtn,
1978            gdt,
1979            mi,
1980            mo,
1981            km,
1982            ibi,
1983            li,
1984            a,
1985            lats=np.array(lats, dtype=np.float32),
1986            lons=np.array(lons, dtype=np.float32),
1987        )
1988        out = go.reshape(newshp)
1989    
1990    elif isinstance(a,tuple):
1991        # Vector
1992        km = a[0].shape[0]
1993        ibi = np.zeros((km), dtype=np.int32)
1994        li = np.zeros(a[0].shape,dtype=np.uint8)
1995        no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret = iplib.interpolate_vector(
1996            method,
1997            method_options,
1998            grid_def_in.gdtn,
1999            np.array(grid_def_in.gdt, dtype=np.int32),
2000            gdtn,
2001            gdt,
2002            mi,
2003            mo,
2004            km,
2005            ibi,
2006            li,
2007            a[0],
2008            a[1],
2009            lats=np.array(lats, dtype=np.float32),
2010            lons=np.array(lons, dtype=np.float32),
2011        )
2012        out = (uo.reshape(newshp),vo.reshape(newshp))
2013
2014    try:
2015        iplib.openmp_set_num_threads(prev_num_threads)
2016    except(AttributeError):
2017        pass
2018
2019    return out