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

29def show_config():
30    """Print grib2io build configuration information."""
31    print(f'grib2io version {__version__} Configuration:\n')
32    print(f'\tg2c library version: {__g2clib_version__}')
33    print(f'\tJPEG compression support: {has_jpeg_support}')
34    print(f'\tPNG compression support: {has_png_support}')
35    print(f'\tAEC compression support: {has_aec_support}')
36    print(f'')
37    print(f'\tNCEP GRIB2 Table Version: {_ncep_grib2_table_version}')

1490def interpolate(a, method: Union[int, str], grid_def_in, grid_def_out, 
1491                method_options=None, num_threads=1):
1492    """
1493    This is the module-level interpolation function.
1494
1495    This interfaces with the grib2io_interp component package that interfaces to
1496    the [NCEPLIBS-ip library](https://github.com/NOAA-EMC/NCEPLIBS-ip).
1497
1498    Parameters
1499    ----------
1500    a : numpy.ndarray or tuple
1501        Input data.  If `a` is a `numpy.ndarray`, scalar interpolation will be
1502        performed.  If `a` is a `tuple`, then vector interpolation will be
1503        performed with the assumption that u = a[0] and v = a[1] and are both
1504        `numpy.ndarray`.
1505
1506        These data are expected to be in 2-dimensional form with shape (ny, nx)
1507        or 3-dimensional (:, ny, nx) where the 1st dimension represents another
1508        spatial, temporal, or classification (i.e. ensemble members) dimension.
1509        The function will properly flatten the (ny,nx) dimensions into (nx * ny)
1510        acceptable for input into the interpolation subroutines.
1511    method
1512        Interpolate method to use. This can either be an integer or string using
1513        the following mapping:
1514
1515        | Interpolate Scheme | Integer Value |
1516        | :---:              | :---:         |
1517        | 'bilinear'         | 0             |
1518        | 'bicubic'          | 1             |
1519        | 'neighbor'         | 2             |
1520        | 'budget'           | 3             |
1521        | 'spectral'         | 4             |
1522        | 'neighbor-budget'  | 6             |
1523
1524    grid_def_in : grib2io.Grib2GridDef
1525        Grib2GridDef object for the input grid.
1526    grid_def_out : grib2io.Grib2GridDef
1527        Grib2GridDef object for the output grid or station points.
1528    method_options : list of ints, optional
1529        Interpolation options. See the NCEPLIBS-ip documentation for
1530        more information on how these are used.
1531    num_threads : int, optional
1532        Number of OpenMP threads to use for interpolation. The default
1533        value is 1. If grib2io_interp was not built with OpenMP, then
1534        this keyword argument and value will have no impact.
1535
1536    Returns
1537    -------
1538    interpolate
1539        Returns a `numpy.ndarray` when scalar interpolation is performed or a
1540        `tuple` of `numpy.ndarray`s when vector interpolation is performed with
1541        the assumptions that 0-index is the interpolated u and 1-index is the
1542        interpolated v.
1543    """
1544    import grib2io_interp
1545    from grib2io_interp import interpolate
1546
1547    prev_num_threads = 1
1548    try:
1549        import grib2io_interp
1550        if grib2io_interp.has_openmp_support:
1551            prev_num_threads = grib2io_interp.get_openmp_threads()
1552            grib2io_interp.set_openmp_threads(num_threads)
1553    except(AttributeError):
1554        pass
1555
1556    if isinstance(method,int) and method not in _interp_schemes.values():
1557        raise ValueError('Invalid interpolation method.')
1558    elif isinstance(method,str):
1559        if method in _interp_schemes.keys():
1560            method = _interp_schemes[method]
1561        else:
1562            raise ValueError('Invalid interpolation method.')
1563
1564    if method_options is None:
1565        method_options = np.zeros((20),dtype=np.int32)
1566        if method in {3,6}:
1567            method_options[0:2] = -1
1568
1569    ni = grid_def_in.npoints
1570    no = grid_def_out.npoints
1571
1572    # Adjust shape of input array(s)
1573    a,newshp = _adjust_array_shape_for_interp(a,grid_def_in,grid_def_out)
1574
1575    # Set lats and lons if stations, else create array for grids.
1576    if grid_def_out.gdtn == -1:
1577        rlat = np.array(grid_def_out.lats,dtype=np.float32)
1578        rlon = np.array(grid_def_out.lons,dtype=np.float32)
1579    else:
1580        rlat = np.zeros((no),dtype=np.float32)
1581        rlon = np.zeros((no),dtype=np.float32)
1582
1583    # Call interpolation subroutines according to type of a.
1584    if isinstance(a,np.ndarray):
1585        # Scalar
1586        if np.any(np.isnan(a)):
1587            ibi = np.zeros((a.shape[0]),dtype=np.int32)+1
1588            li = np.where(np.isnan(a),0,1).astype(np.int8)
1589        else:
1590            ibi = np.zeros((a.shape[0]),dtype=np.int32)
1591            li = np.zeros(a.shape,dtype=np.int8)
1592        go = np.zeros((a.shape[0],no),dtype=np.float32)
1593        no,ibo,lo,iret = interpolate.interpolate_scalar(method,method_options,
1594                                                 grid_def_in.gdtn,grid_def_in.gdt,
1595                                                 grid_def_out.gdtn,grid_def_out.gdt,
1596                                                 ibi,li.T,a.T,go.T,rlat,rlon)
1597        lo = lo.T.reshape(newshp)
1598        out = go.reshape(newshp)
1599        out = np.where(lo==0,np.nan,out)
1600    elif isinstance(a,tuple):
1601        # Vector
1602        if np.any(np.isnan(a)):
1603            ibi = np.zeros((a.shape[0]),dtype=np.int32)+1
1604            li = np.where(np.isnan(a),0,1).astype(np.int8)
1605        else:
1606            ibi = np.zeros((a.shape[0]),dtype=np.int32)
1607            li = np.zeros(a.shape,dtype=np.int8)
1608        uo = np.zeros((a[0].shape[0],no),dtype=np.float32)
1609        vo = np.zeros((a[1].shape[0],no),dtype=np.float32)
1610        crot = np.ones((no),dtype=np.float32)
1611        srot = np.zeros((no),dtype=np.float32)
1612        no,ibo,lo,iret = interpolate.interpolate_vector(method,method_options,
1613                                                 grid_def_in.gdtn,grid_def_in.gdt,
1614                                                 grid_def_out.gdtn,grid_def_out.gdt,
1615                                                 ibi,li.T,a[0].T,a[1].T,uo.T,vo.T,
1616                                                 rlat,rlon,crot,srot)
1617        del crot
1618        del srot
1619        lo = lo[:,0].reshape(newshp)
1620        uo = uo.reshape(new)
1621        vo = vo.reshape(new)
1622        uo = np.where(lo==0,np.nan,uo)
1623        vo = np.where(lo==0,np.nan,vo)
1624        out = (uo,vo)
1625
1626    del rlat
1627    del rlon
1628
1629    try:
1630        if grib2io_interp.has_openmp_support:
1631            grib2io_interp.set_openmp_threads(prev_num_threads)
1632    except(AttributeError):
1633        pass
1634
1635    return out

1638def interpolate_to_stations(a, method, grid_def_in, lats, lons,
1639                            method_options=None, num_threads=1):
1640    """
1641    Module-level interpolation function for interpolation to stations.
1642
1643    Interfaces with the grib2io_interp component package that interfaces to the
1644    [NCEPLIBS-ip
1645    library](https://github.com/NOAA-EMC/NCEPLIBS-ip). It supports scalar and
1646    vector interpolation according to the type of object a.
1647
1648    Parameters
1649    ----------
1650    a : numpy.ndarray or tuple
1651        Input data.  If `a` is a `numpy.ndarray`, scalar interpolation will be
1652        performed.  If `a` is a `tuple`, then vector interpolation will be
1653        performed with the assumption that u = a[0] and v = a[1] and are both
1654        `numpy.ndarray`.
1655
1656        These data are expected to be in 2-dimensional form with shape (ny, nx)
1657        or 3-dimensional (:, ny, nx) where the 1st dimension represents another
1658        spatial, temporal, or classification (i.e. ensemble members) dimension.
1659        The function will properly flatten the (ny,nx) dimensions into (nx * ny)
1660        acceptable for input into the interpolation subroutines.
1661    method
1662        Interpolate method to use. This can either be an integer or string using
1663        the following mapping:
1664
1665        | Interpolate Scheme | Integer Value |
1666        | :---:              | :---:         |
1667        | 'bilinear'         | 0             |
1668        | 'bicubic'          | 1             |
1669        | 'neighbor'         | 2             |
1670        | 'budget'           | 3             |
1671        | 'spectral'         | 4             |
1672        | 'neighbor-budget'  | 6             |
1673
1674    grid_def_in : grib2io.Grib2GridDef
1675        Grib2GridDef object for the input grid.
1676    lats : numpy.ndarray or list
1677        Latitudes for station points
1678    lons : numpy.ndarray or list
1679        Longitudes for station points
1680    method_options : list of ints, optional
1681        Interpolation options. See the NCEPLIBS-ip documentation for
1682        more information on how these are used.
1683    num_threads : int, optional
1684        Number of OpenMP threads to use for interpolation. The default
1685        value is 1. If grib2io_interp was not built with OpenMP, then
1686        this keyword argument and value will have no impact.
1687
1688    Returns
1689    -------
1690    interpolate_to_stations
1691        Returns a `numpy.ndarray` when scalar interpolation is performed or a
1692        `tuple` of `numpy.ndarray`s when vector interpolation is performed with
1693        the assumptions that 0-index is the interpolated u and 1-index is the
1694        interpolated v.
1695    """
1696    import grib2io_interp
1697    from grib2io_interp import interpolate
1698
1699    prev_num_threads = 1
1700    try:
1701        import grib2io_interp
1702        if grib2io_interp.has_openmp_support:
1703            prev_num_threads = grib2io_interp.get_openmp_threads()
1704            grib2io_interp.set_openmp_threads(num_threads)
1705    except(AttributeError):
1706        pass
1707
1708    if isinstance(method,int) and method not in _interp_schemes.values():
1709        raise ValueError('Invalid interpolation method.')
1710    elif isinstance(method,str):
1711        if method in _interp_schemes.keys():
1712            method = _interp_schemes[method]
1713        else:
1714            raise ValueError('Invalid interpolation method.')
1715
1716    if method_options is None:
1717        method_options = np.zeros((20),dtype=np.int32)
1718        if method in {3,6}:
1719            method_options[0:2] = -1
1720
1721    # Check lats and lons
1722    if isinstance(lats,list):
1723        nlats = len(lats)
1724    elif isinstance(lats,np.ndarray) and len(lats.shape) == 1:
1725        nlats = lats.shape[0]
1726    else:
1727        raise ValueError('Station latitudes must be a list or 1-D NumPy array.')
1728    if isinstance(lons,list):
1729        nlons = len(lons)
1730    elif isinstance(lons,np.ndarray) and len(lons.shape) == 1:
1731        nlons = lons.shape[0]
1732    else:
1733        raise ValueError('Station longitudes must be a list or 1-D NumPy array.')
1734    if nlats != nlons:
1735        raise ValueError('Station lats and lons must be same size.')
1736
1737    ni = grid_def_in.npoints
1738    no = nlats
1739
1740    # Adjust shape of input array(s)
1741    a,newshp = _adjust_array_shape_for_interp_stations(a,grid_def_in,no)
1742
1743    # Set lats and lons if stations
1744    rlat = np.array(lats,dtype=np.float32)
1745    rlon = np.array(lons,dtype=np.float32)
1746
1747    # Use gdtn = -1 for stations and an empty template array
1748    gdtn = -1
1749    gdt = np.zeros((200),dtype=np.int32)
1750
1751    # Call interpolation subroutines according to type of a.
1752    if isinstance(a,np.ndarray):
1753        # Scalar
1754        ibi = np.zeros((a.shape[0]),dtype=np.int32)
1755        li = np.zeros(a.shape,dtype=np.int32)
1756        go = np.zeros((a.shape[0],no),dtype=np.float32)
1757        no,ibo,lo,iret = interpolate.interpolate_scalar(method,method_options,
1758                                                 grid_def_in.gdtn,grid_def_in.gdt,
1759                                                 gdtn,gdt,
1760                                                 ibi,li.T,a.T,go.T,rlat,rlon)
1761        out = go.reshape(newshp)
1762    elif isinstance(a,tuple):
1763        # Vector
1764        ibi = np.zeros((a[0].shape[0]),dtype=np.int32)
1765        li = np.zeros(a[0].shape,dtype=np.int32)
1766        uo = np.zeros((a[0].shape[0],no),dtype=np.float32)
1767        vo = np.zeros((a[1].shape[0],no),dtype=np.float32)
1768        crot = np.ones((no),dtype=np.float32)
1769        srot = np.zeros((no),dtype=np.float32)
1770        no,ibo,lo,iret = interpolate.interpolate_vector(method,method_options,
1771                                                 grid_def_in.gdtn,grid_def_in.gdt,
1772                                                 gdtn,gdt,
1773                                                 ibi,li.T,a[0].T,a[1].T,uo.T,vo.T,
1774                                                 rlat,rlon,crot,srot)
1775        del crot
1776        del srot
1777        out = (uo.reshape(newshp),vo.reshape(newshp))
1778
1779    del rlat
1780    del rlon
1781
1782    try:
1783        if grib2io_interp.has_openmp_support:
1784            grib2io_interp.set_openmp_threads(prev_num_threads)
1785    except(AttributeError):
1786        pass
1787
1788    return out