grib2io
Introduction
grib2io is a Python package that provides an interface to the NCEP GRIB2 C (g2c) library for the purpose of reading and writing WMO GRIdded Binary, Edition 2 (GRIB2) messages. A physical file can contain one or more GRIB2 messages.
GRIB2 file IO is performed directly in Python. The unpacking/packing of GRIB2 integer, coded metadata and data sections is performed by the g2c library functions via the g2clib Cython extension module. The decoding/encoding of GRIB2 metadata is translated into more descriptive, plain language metadata by providing a mapping of the integer coded metadata to the appropriate GRIB2 code tables. The NCEP GRIB2 code tables are a part of the grib2io package.
Interpolation
As of grib2io v2.4.0, spatial interpolation via NCEPLIPS-ip Fortran library is now a part of the grib2io package. The separate component package, grib2io-interp, has been deprecated. grib2io-interp provided interpolation via F2PY interface to NCEPLIBS-ip, which has become difficult since the removal of distutils from Python 3.12+.
NCEPLIBS-ip interpolation Fortran subroutines contain the BIND(C)
attribute which
provides an equivalent C-interface. grib2io now provides a Cython-based interface, iplib
,
to these Fortran subroutines via their C-interface. If NCEPLIBS-ip was built with
OpenMP support, iplib
will provide functions for getting and setting the number of
OpenMP threads.
Xarray Backend
grib2io provides a Xarray backend engine so that many GRIB2 messages can be represented as N-dimensional DataArray objects and collected along common coordinates as Datasets or DataTrees. The Xarray backend engine API is experimental and is subject to change without backward compatibility.
Tutorials
The following Jupyter Notebooks are available as tutorials:
1from ._grib2io import * 2from ._grib2io import __doc__ 3from ._grib2io import _Grib2Message 4 5__all__ = ['open', 'show_config', 'interpolate', 'interpolate_to_stations', 6 'tables', 'templates', 'utils', 7 'Grib2Message', '_Grib2Message', 'Grib2GridDef', 'msgs_from_index'] 8 9try: 10 from . import __config__ 11 __version__ = __config__.grib2io_version 12 has_interpolation = __config__.has_interpolation 13 has_openmp_support = __config__.has_openmp_support 14 g2c_static = __config__.g2c_static 15 ip_static = __config__.ip_static 16 extra_objects = __config__.extra_objects 17except(ImportError): 18 pass 19 20from .g2clib import __version__ as __g2clib_version__ 21from .g2clib import _has_jpeg 22from .g2clib import _has_png 23from .g2clib import _has_aec 24 25has_jpeg_support = bool(_has_jpeg) 26has_png_support = bool(_has_png) 27has_aec_support = bool(_has_aec) 28 29from .tables.originating_centers import _ncep_grib2_table_version 30ncep_grib2_table_version = _ncep_grib2_table_version 31g2c_version = __g2clib_version__ 32 33def show_config(): 34 """Print grib2io build configuration information.""" 35 print(f'grib2io version {__version__} Configuration:') 36 print(f'') 37 print(f'NCEPLIBS-g2c library version: {__g2clib_version__}') 38 print(f'\tStatic library: {g2c_static}') 39 print(f'\tJPEG compression support: {has_jpeg_support}') 40 print(f'\tPNG compression support: {has_png_support}') 41 print(f'\tAEC compression support: {has_aec_support}') 42 print(f'') 43 print(f'NCEPLIPS-ip support: {has_interpolation}') 44 print(f'\tStatic library: {ip_static}') 45 print(f'\tOpenMP support: {has_openmp_support}') 46 print(f'') 47 print(f'Static libs:') 48 for lib in extra_objects: 49 print(f'\t{lib}') 50 print(f'') 51 print(f'NCEP GRIB2 Table Version: {_ncep_grib2_table_version}')
95class open(): 96 """ 97 GRIB2 File Object. 98 99 This class can accommodate a physical file with one or more GRIB2 100 messages or a bytes object containing a GRIB2 messages. 101 102 A physical file can contain one or more GRIB2 messages. When instantiated, 103 class `grib2io.open`, the file named `filename` is opened for reading (`mode 104 = 'r'`) and is automatically indexed. The indexing procedure reads some of 105 the GRIB2 metadata for all GRIB2 Messages. A GRIB2 Message may contain 106 submessages whereby Section 2-7 can be repeated. grib2io accommodates for 107 this by flattening any GRIB2 submessages into multiple individual messages. 108 109 It is important to note that GRIB2 files from some Meteorological agencies 110 contain other data than GRIB2 messages. GRIB2 files from ECMWF can contain 111 GRIB1 and GRIB2 messages. grib2io checks for these and safely ignores them. 112 113 Attributes 114 ---------- 115 closed : bool 116 `True` is file handle is close; `False` otherwise. 117 current_message : int 118 Current position of the file in units of GRIB2 Messages. (read only) 119 levels : tuple 120 Tuple containing a unique list of wgrib2-formatted level/layer strings. 121 messages : int 122 Count of GRIB2 Messages contained in the file. 123 mode : str 124 File IO mode of opening the file. 125 name : str 126 Full path name of the GRIB2 file. 127 save_index : bool 128 Whether to save a pickle-based index file for the GRIB2 file. Default is `True`. 129 size : int 130 Size of the file in units of bytes. 131 use_index 132 Whether to use an existing pickle-based index file for the GRIB2 file. Default is `True`. 133 variables : tuple 134 Tuple containing a unique list of variable short names (i.e. GRIB2 135 abbreviation names). 136 """ 137 138 __slots__ = ('_fileid', '_filehandle', '_hasindex', '_index', 139 '_pos', 'closed', 'current_message', 'messages', 'mode', 140 'name', 'size', 'save_index', 'use_index', '_msgs') 141 142 def __init__( 143 self, 144 filename: Union[bytes, str, Path], 145 mode: Literal["r", "w", "x"] = "r", 146 *, 147 save_index = True, 148 use_index = True, 149 _xarray_backend = False, 150 **kwargs, 151 ): 152 """ 153 Initialize GRIB2 File object instance. 154 155 Parameters 156 ---------- 157 filename 158 File path containing GRIB2 messages OR bytes. 159 mode 160 File access mode where "r" opens the files for reading only; "w" 161 opens the file for overwriting and "x" for writing to a new file. 162 save_index 163 Whether to save a pickle-based index file for the GRIB2 file. Default is True. 164 use_index 165 Whether to use an existing pickle-based index file for the GRIB2 file. Default is True. 166 """ 167 168 # All write modes are read/write. 169 # All modes are binary. 170 if mode in ("a", "x", "w"): 171 mode += "+" 172 mode = mode + "b" 173 174 self.mode = mode 175 self.save_index = save_index 176 self.use_index = use_index 177 178 if isinstance(filename, bytes): 179 if 'r' not in self.mode: 180 raise ValueError("bytes grib2 can only be opened as read") 181 self.current_message = 0 182 if filename[:2] == _GZIP_HEADER: 183 filename = gzip.decompress(filename) 184 if filename[:4]+filename[-4:] != b'GRIB7777': 185 raise ValueError("Invalid GRIB bytes") 186 self._filehandle = BytesIO(filename) 187 self.name = "<in-memory-file>" 188 self.size = len(filename) 189 self._fileid = hashlib.sha1((self.name+str(self.size)).encode('ASCII')).hexdigest() 190 self._index = build_index(self._filehandle) 191 self._msgs = msgs_from_index(self._index, filehandle=self._filehandle) 192 self.messages = len(self._msgs) 193 194 else: 195 self.current_message = 0 196 if 'r' in mode: 197 self._filehandle = builtins.open(filename, mode=mode) 198 # Some GRIB2 files are gzipped, so check for that here, but 199 # raise error when using xarray backend. 200 # Gzip files contain a 2-byte header b'\x1f\x8b'. 201 if self._filehandle.read(2) == _GZIP_HEADER: 202 self._filehandle.close() 203 if _xarray_backend: 204 raise RuntimeError('Gzip GRIB2 files are not supported by the Xarray backend.') 205 self._filehandle = gzip.open(filename, mode=mode) 206 else: 207 self._filehandle = builtins.open(filename, mode=mode) 208 else: 209 self._filehandle = builtins.open(filename, mode=mode) 210 self.name = os.path.abspath(filename) 211 self.name = os.path.abspath(filename) 212 fstat = os.stat(self.name) 213 self.size = fstat.st_size 214 self._fileid = hashlib.sha1((self.name+str(fstat.st_ino)+ 215 str(self.size)).encode('ASCII')).hexdigest() 216 217 if 'r' in self.mode: 218 219 indexfile = f"{self.name}_{self._fileid}.grib2ioidx" 220 if self.use_index and os.path.exists(indexfile): 221 try: 222 with builtins.open(indexfile, 'rb') as file: 223 index = pickle.load(file) 224 self._index = index 225 except Exception as e: 226 print(f"found indexfile: {indexfile}, but unable to load it: {e}") 227 print(f"re-forming index from grib2file, but not writing indexfile") 228 self._index = build_index(self._filehandle) 229 else: 230 self._index = build_index(self._filehandle) 231 if self.save_index: 232 try: 233 serialize_index(self._index, indexfile) 234 except Exception as e: 235 print(f"index was not serialized for future use: {e}") 236 237 self._msgs = msgs_from_index(self._index, filehandle=self._filehandle) 238 239 self.messages = len(self._msgs) 240 elif 'w' or 'x' in self.mode: 241 self.messages = 0 242 self.current_message = None 243 244 self.closed = self._filehandle.closed 245 246 247 def __delete__(self, instance): 248 self.close() 249 del self._index 250 251 252 def __enter__(self): 253 return self 254 255 256 def __exit__(self, atype, value, traceback): 257 self.close() 258 259 260 def __iter__(self): 261 yield from self._msgs 262 263 264 def __len__(self): 265 return self.messages 266 267 268 def __repr__(self): 269 strings = [] 270 for k in self.__slots__: 271 if k.startswith('_'): continue 272 strings.append('%s = %s\n'%(k,eval('self.'+k))) 273 return ''.join(strings) 274 275 276 def __getitem__(self, key): 277 if isinstance(key,int): 278 if abs(key) >= len(self._msgs): 279 raise IndexError("index out of range") 280 else: 281 return self._msgs[key] 282 elif isinstance(key,str): 283 return self.select(shortName=key) 284 elif isinstance(key,slice): 285 return self._msgs[key] 286 else: 287 raise KeyError('Key must be an integer, slice, or GRIB2 variable shortName.') 288 289 290 291 @property 292 def levels(self): 293 """Provides a unique tuple of level strings.""" 294 if self._hasindex: 295 return tuple(sorted(set([msg.level for msg in self._msgs]))) 296 else: 297 return None 298 299 300 @property 301 def variables(self): 302 """Provides a unique tuple of variable shortName strings.""" 303 if self._hasindex: 304 return tuple(sorted(set([msg.shortName for msg in self._msgs]))) 305 else: 306 return None 307 308 309 def close(self): 310 """Close the file handle.""" 311 if not self._filehandle.closed: 312 self.messages = 0 313 self.current_message = 0 314 self._filehandle.close() 315 self.closed = self._filehandle.closed 316 317 318 def read(self, size: Optional[int]=None): 319 """ 320 Read size amount of GRIB2 messages from the current position. 321 322 If no argument is given, then size is None and all messages are returned 323 from the current position in the file. This read method follows the 324 behavior of Python's builtin open() function, but whereas that operates 325 on units of bytes, we operate on units of GRIB2 messages. 326 327 Parameters 328 ---------- 329 size: default=None 330 The number of GRIB2 messages to read from the current position. If 331 no argument is give, the default value is None and remainder of 332 the file is read. 333 334 Returns 335 ------- 336 read 337 ``Grib2Message`` object when size = 1 or a list of Grib2Messages 338 when size > 1. 339 """ 340 if size is not None and size < 0: 341 size = None 342 if size is None or size > 1: 343 start = self.tell() 344 stop = self.messages if size is None else start+size 345 if size is None: 346 self.current_message = self.messages-1 347 else: 348 self.current_message += size 349 return self._msgs[slice(start,stop,1)] 350 elif size == 1: 351 self.current_message += 1 352 return self._msgs[self.current_message] 353 else: 354 None 355 356 357 def seek(self, pos: int): 358 """ 359 Set the position within the file in units of GRIB2 messages. 360 361 Parameters 362 ---------- 363 pos 364 The GRIB2 Message number to set the file pointer to. 365 """ 366 if self._hasindex: 367 self._filehandle.seek(self._index['sectionOffset'][0][pos]) 368 self.current_message = pos 369 370 371 def tell(self): 372 """Returns the position of the file in units of GRIB2 Messages.""" 373 return self.current_message 374 375 376 def select(self, **kwargs): 377 """Select GRIB2 messages by `Grib2Message` attributes.""" 378 # TODO: Added ability to process multiple values for each keyword (attribute) 379 idxs = [] 380 nkeys = len(kwargs.keys()) 381 for k,v in kwargs.items(): 382 for m in self._msgs: 383 if hasattr(m,k) and getattr(m,k) == v: idxs.append(m._msgnum) 384 idxs = np.array(idxs,dtype=np.int32) 385 return [self._msgs[i] for i in [ii[0] for ii in collections.Counter(idxs).most_common() if ii[1] == nkeys]] 386 387 388 def write(self, msg): 389 """ 390 Writes GRIB2 message object to file. 391 392 Parameters 393 ---------- 394 msg 395 GRIB2 message objects to write to file. 396 """ 397 if isinstance(msg,list): 398 for m in msg: 399 self.write(m) 400 return 401 402 if issubclass(msg.__class__,_Grib2Message): 403 # TODO: We can consider letting pack return packed bytes instead of associating with message object 404 if hasattr(msg,'_msg'): 405 # write already packed bytes 406 self._filehandle.write(msg._msg) 407 else: 408 if msg._signature == msg._generate_signature() and \ 409 msg._data is None and \ 410 hasattr(msg._ondiskarray,'filehandle'): 411 # write unchanged message from input 412 offset = msg._ondiskarray.filehandle.tell() 413 msg._ondiskarray.filehandle.seek(msg._ondiskarray.offset) 414 self._filehandle.write(msg._ondiskarray.filehandle.read(msg.section0[-1])) 415 msg._ondiskarray.filehandle.seek(offset) 416 else: 417 msg.pack() 418 self._filehandle.write(msg._msg) 419 self.size = os.path.getsize(self.name) 420 self.messages += 1 421 else: 422 raise TypeError("msg must be a Grib2Message object.") 423 return 424 425 426 def flush(self): 427 """Flush the file object buffer.""" 428 self._filehandle.flush() 429 430 431 def levels_by_var(self, name: str): 432 """ 433 Return a list of level strings given a variable shortName. 434 435 Parameters 436 ---------- 437 name 438 Grib2Message variable shortName 439 440 Returns 441 ------- 442 levels_by_var 443 A list of unique level strings. 444 """ 445 return list(sorted(set([msg.level for msg in self.select(shortName=name)]))) 446 447 448 def vars_by_level(self, level: str): 449 """ 450 Return a list of variable shortName strings given a level. 451 452 Parameters 453 ---------- 454 level 455 Grib2Message variable level 456 457 Returns 458 ------- 459 vars_by_level 460 A list of unique variable shortName strings. 461 """ 462 return list(sorted(set([msg.shortName for msg in self.select(level=level)])))
GRIB2 File Object.
This class can accommodate a physical file with one or more GRIB2 messages or a bytes object containing a GRIB2 messages.
A physical file can contain one or more GRIB2 messages. When instantiated,
class grib2io.open
, the file named filename
is opened for reading (mode
= 'r'
) and is automatically indexed. The indexing procedure reads some of
the GRIB2 metadata for all GRIB2 Messages. A GRIB2 Message may contain
submessages whereby Section 2-7 can be repeated. grib2io accommodates for
this by flattening any GRIB2 submessages into multiple individual messages.
It is important to note that GRIB2 files from some Meteorological agencies contain other data than GRIB2 messages. GRIB2 files from ECMWF can contain GRIB1 and GRIB2 messages. grib2io checks for these and safely ignores them.
Attributes
- closed (bool):
True
is file handle is close;False
otherwise. - current_message (int): Current position of the file in units of GRIB2 Messages. (read only)
- levels (tuple): Tuple containing a unique list of wgrib2-formatted level/layer strings.
- messages (int): Count of GRIB2 Messages contained in the file.
- mode (str): File IO mode of opening the file.
- name (str): Full path name of the GRIB2 file.
- save_index (bool):
Whether to save a pickle-based index file for the GRIB2 file. Default is
True
. - size (int): Size of the file in units of bytes.
- use_index: Whether to use an existing pickle-based index file for the GRIB2 file. Default is
True
. - variables (tuple): Tuple containing a unique list of variable short names (i.e. GRIB2 abbreviation names).
142 def __init__( 143 self, 144 filename: Union[bytes, str, Path], 145 mode: Literal["r", "w", "x"] = "r", 146 *, 147 save_index = True, 148 use_index = True, 149 _xarray_backend = False, 150 **kwargs, 151 ): 152 """ 153 Initialize GRIB2 File object instance. 154 155 Parameters 156 ---------- 157 filename 158 File path containing GRIB2 messages OR bytes. 159 mode 160 File access mode where "r" opens the files for reading only; "w" 161 opens the file for overwriting and "x" for writing to a new file. 162 save_index 163 Whether to save a pickle-based index file for the GRIB2 file. Default is True. 164 use_index 165 Whether to use an existing pickle-based index file for the GRIB2 file. Default is True. 166 """ 167 168 # All write modes are read/write. 169 # All modes are binary. 170 if mode in ("a", "x", "w"): 171 mode += "+" 172 mode = mode + "b" 173 174 self.mode = mode 175 self.save_index = save_index 176 self.use_index = use_index 177 178 if isinstance(filename, bytes): 179 if 'r' not in self.mode: 180 raise ValueError("bytes grib2 can only be opened as read") 181 self.current_message = 0 182 if filename[:2] == _GZIP_HEADER: 183 filename = gzip.decompress(filename) 184 if filename[:4]+filename[-4:] != b'GRIB7777': 185 raise ValueError("Invalid GRIB bytes") 186 self._filehandle = BytesIO(filename) 187 self.name = "<in-memory-file>" 188 self.size = len(filename) 189 self._fileid = hashlib.sha1((self.name+str(self.size)).encode('ASCII')).hexdigest() 190 self._index = build_index(self._filehandle) 191 self._msgs = msgs_from_index(self._index, filehandle=self._filehandle) 192 self.messages = len(self._msgs) 193 194 else: 195 self.current_message = 0 196 if 'r' in mode: 197 self._filehandle = builtins.open(filename, mode=mode) 198 # Some GRIB2 files are gzipped, so check for that here, but 199 # raise error when using xarray backend. 200 # Gzip files contain a 2-byte header b'\x1f\x8b'. 201 if self._filehandle.read(2) == _GZIP_HEADER: 202 self._filehandle.close() 203 if _xarray_backend: 204 raise RuntimeError('Gzip GRIB2 files are not supported by the Xarray backend.') 205 self._filehandle = gzip.open(filename, mode=mode) 206 else: 207 self._filehandle = builtins.open(filename, mode=mode) 208 else: 209 self._filehandle = builtins.open(filename, mode=mode) 210 self.name = os.path.abspath(filename) 211 self.name = os.path.abspath(filename) 212 fstat = os.stat(self.name) 213 self.size = fstat.st_size 214 self._fileid = hashlib.sha1((self.name+str(fstat.st_ino)+ 215 str(self.size)).encode('ASCII')).hexdigest() 216 217 if 'r' in self.mode: 218 219 indexfile = f"{self.name}_{self._fileid}.grib2ioidx" 220 if self.use_index and os.path.exists(indexfile): 221 try: 222 with builtins.open(indexfile, 'rb') as file: 223 index = pickle.load(file) 224 self._index = index 225 except Exception as e: 226 print(f"found indexfile: {indexfile}, but unable to load it: {e}") 227 print(f"re-forming index from grib2file, but not writing indexfile") 228 self._index = build_index(self._filehandle) 229 else: 230 self._index = build_index(self._filehandle) 231 if self.save_index: 232 try: 233 serialize_index(self._index, indexfile) 234 except Exception as e: 235 print(f"index was not serialized for future use: {e}") 236 237 self._msgs = msgs_from_index(self._index, filehandle=self._filehandle) 238 239 self.messages = len(self._msgs) 240 elif 'w' or 'x' in self.mode: 241 self.messages = 0 242 self.current_message = None 243 244 self.closed = self._filehandle.closed
Initialize GRIB2 File object instance.
Parameters
- filename: File path containing GRIB2 messages OR bytes.
- mode: File access mode where "r" opens the files for reading only; "w" opens the file for overwriting and "x" for writing to a new file.
- save_index: Whether to save a pickle-based index file for the GRIB2 file. Default is True.
- use_index: Whether to use an existing pickle-based index file for the GRIB2 file. Default is True.
291 @property 292 def levels(self): 293 """Provides a unique tuple of level strings.""" 294 if self._hasindex: 295 return tuple(sorted(set([msg.level for msg in self._msgs]))) 296 else: 297 return None
Provides a unique tuple of level strings.
300 @property 301 def variables(self): 302 """Provides a unique tuple of variable shortName strings.""" 303 if self._hasindex: 304 return tuple(sorted(set([msg.shortName for msg in self._msgs]))) 305 else: 306 return None
Provides a unique tuple of variable shortName strings.
309 def close(self): 310 """Close the file handle.""" 311 if not self._filehandle.closed: 312 self.messages = 0 313 self.current_message = 0 314 self._filehandle.close() 315 self.closed = self._filehandle.closed
Close the file handle.
318 def read(self, size: Optional[int]=None): 319 """ 320 Read size amount of GRIB2 messages from the current position. 321 322 If no argument is given, then size is None and all messages are returned 323 from the current position in the file. This read method follows the 324 behavior of Python's builtin open() function, but whereas that operates 325 on units of bytes, we operate on units of GRIB2 messages. 326 327 Parameters 328 ---------- 329 size: default=None 330 The number of GRIB2 messages to read from the current position. If 331 no argument is give, the default value is None and remainder of 332 the file is read. 333 334 Returns 335 ------- 336 read 337 ``Grib2Message`` object when size = 1 or a list of Grib2Messages 338 when size > 1. 339 """ 340 if size is not None and size < 0: 341 size = None 342 if size is None or size > 1: 343 start = self.tell() 344 stop = self.messages if size is None else start+size 345 if size is None: 346 self.current_message = self.messages-1 347 else: 348 self.current_message += size 349 return self._msgs[slice(start,stop,1)] 350 elif size == 1: 351 self.current_message += 1 352 return self._msgs[self.current_message] 353 else: 354 None
Read size amount of GRIB2 messages from the current position.
If no argument is given, then size is None and all messages are returned from the current position in the file. This read method follows the behavior of Python's builtin open() function, but whereas that operates on units of bytes, we operate on units of GRIB2 messages.
Parameters
- size (default=None): The number of GRIB2 messages to read from the current position. If no argument is give, the default value is None and remainder of the file is read.
Returns
- read:
Grib2Message
object when size = 1 or a list of Grib2Messages when size > 1.
357 def seek(self, pos: int): 358 """ 359 Set the position within the file in units of GRIB2 messages. 360 361 Parameters 362 ---------- 363 pos 364 The GRIB2 Message number to set the file pointer to. 365 """ 366 if self._hasindex: 367 self._filehandle.seek(self._index['sectionOffset'][0][pos]) 368 self.current_message = pos
Set the position within the file in units of GRIB2 messages.
Parameters
- pos: The GRIB2 Message number to set the file pointer to.
371 def tell(self): 372 """Returns the position of the file in units of GRIB2 Messages.""" 373 return self.current_message
Returns the position of the file in units of GRIB2 Messages.
376 def select(self, **kwargs): 377 """Select GRIB2 messages by `Grib2Message` attributes.""" 378 # TODO: Added ability to process multiple values for each keyword (attribute) 379 idxs = [] 380 nkeys = len(kwargs.keys()) 381 for k,v in kwargs.items(): 382 for m in self._msgs: 383 if hasattr(m,k) and getattr(m,k) == v: idxs.append(m._msgnum) 384 idxs = np.array(idxs,dtype=np.int32) 385 return [self._msgs[i] for i in [ii[0] for ii in collections.Counter(idxs).most_common() if ii[1] == nkeys]]
Select GRIB2 messages by Grib2Message
attributes.
388 def write(self, msg): 389 """ 390 Writes GRIB2 message object to file. 391 392 Parameters 393 ---------- 394 msg 395 GRIB2 message objects to write to file. 396 """ 397 if isinstance(msg,list): 398 for m in msg: 399 self.write(m) 400 return 401 402 if issubclass(msg.__class__,_Grib2Message): 403 # TODO: We can consider letting pack return packed bytes instead of associating with message object 404 if hasattr(msg,'_msg'): 405 # write already packed bytes 406 self._filehandle.write(msg._msg) 407 else: 408 if msg._signature == msg._generate_signature() and \ 409 msg._data is None and \ 410 hasattr(msg._ondiskarray,'filehandle'): 411 # write unchanged message from input 412 offset = msg._ondiskarray.filehandle.tell() 413 msg._ondiskarray.filehandle.seek(msg._ondiskarray.offset) 414 self._filehandle.write(msg._ondiskarray.filehandle.read(msg.section0[-1])) 415 msg._ondiskarray.filehandle.seek(offset) 416 else: 417 msg.pack() 418 self._filehandle.write(msg._msg) 419 self.size = os.path.getsize(self.name) 420 self.messages += 1 421 else: 422 raise TypeError("msg must be a Grib2Message object.") 423 return
Writes GRIB2 message object to file.
Parameters
- msg: GRIB2 message objects to write to file.
431 def levels_by_var(self, name: str): 432 """ 433 Return a list of level strings given a variable shortName. 434 435 Parameters 436 ---------- 437 name 438 Grib2Message variable shortName 439 440 Returns 441 ------- 442 levels_by_var 443 A list of unique level strings. 444 """ 445 return list(sorted(set([msg.level for msg in self.select(shortName=name)])))
Return a list of level strings given a variable shortName.
Parameters
- name: Grib2Message variable shortName
Returns
- levels_by_var: A list of unique level strings.
448 def vars_by_level(self, level: str): 449 """ 450 Return a list of variable shortName strings given a level. 451 452 Parameters 453 ---------- 454 level 455 Grib2Message variable level 456 457 Returns 458 ------- 459 vars_by_level 460 A list of unique variable shortName strings. 461 """ 462 return list(sorted(set([msg.shortName for msg in self.select(level=level)])))
Return a list of variable shortName strings given a level.
Parameters
- level: Grib2Message variable level
Returns
- vars_by_level: A list of unique variable shortName strings.
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}')
Print grib2io build configuration information.
1795def interpolate(a, 1796 method: Union[int, str], 1797 grid_def_in, 1798 grid_def_out, 1799 method_options=None, 1800 num_threads=1, 1801): 1802 """ 1803 This is the module-level interpolation function. 1804 1805 This interfaces with the 4-byte (32-bit float) [NCEPLIBS-ip library](https://github.com/NOAA-EMC/NCEPLIBS-ip) 1806 through grib2io's internal iplib Cython extension module. 1807 1808 Parameters 1809 ---------- 1810 a : numpy.ndarray or tuple 1811 Input data. If `a` is a `numpy.ndarray`, scalar interpolation will be 1812 performed. If `a` is a `tuple`, then vector interpolation will be 1813 performed with the assumption that u = a[0] and v = a[1] and are both 1814 `numpy.ndarray`. 1815 1816 These data are expected to be in 2-dimensional form with shape (ny, nx) 1817 or 3-dimensional (:, ny, nx) where the 1st dimension represents another 1818 spatial, temporal, or classification (i.e. ensemble members) dimension. 1819 The function will properly flatten the (ny,nx) dimensions into (nx * ny) 1820 acceptable for input into the interpolation subroutines. If needed, these 1821 data will be converted to `np.float32`. 1822 method 1823 Interpolate method to use. This can either be an integer or string using 1824 the following mapping: 1825 1826 | Interpolate Scheme | Integer Value | 1827 | :---: | :---: | 1828 | 'bilinear' | 0 | 1829 | 'bicubic' | 1 | 1830 | 'neighbor' | 2 | 1831 | 'budget' | 3 | 1832 | 'spectral' | 4 | 1833 | 'neighbor-budget' | 6 | 1834 1835 grid_def_in : grib2io.Grib2GridDef 1836 Grib2GridDef object for the input grid. 1837 grid_def_out : grib2io.Grib2GridDef 1838 Grib2GridDef object for the output grid or station points. 1839 method_options : list of ints, optional 1840 Interpolation options. See the NCEPLIBS-ip documentation for 1841 more information on how these are used. 1842 num_threads : int, optional 1843 Number of OpenMP threads to use for interpolation. The default 1844 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 1845 was not built with OpenMP, then this keyword argument and value 1846 will have no impact. 1847 1848 Returns 1849 ------- 1850 interpolate 1851 Returns a `numpy.ndarray` of dtype `np.float32` when scalar interpolation 1852 is performed or a `tuple` of `numpy.ndarray`s when vector interpolation is 1853 performed with the assumptions that 0-index is the interpolated u and 1854 1-index is the interpolated v. 1855 """ 1856 1857 from . import iplib 1858 1859 prev_num_threads = 1 1860 try: 1861 prev_num_threads = iplib.openmp_get_num_threads() 1862 iplib.openmp_set_num_threads(num_threads) 1863 except(AttributeError): 1864 pass 1865 1866 print(f"grib2io.interpolate thread report: OpenMP num threads = {iplib.openmp_get_num_threads()}") 1867 1868 if isinstance(method,int) and method not in _interp_schemes.values(): 1869 raise ValueError('Invalid interpolation method.') 1870 elif isinstance(method,str): 1871 if method in _interp_schemes.keys(): 1872 method = _interp_schemes[method] 1873 else: 1874 raise ValueError('Invalid interpolation method.') 1875 1876 if method_options is None: 1877 method_options = np.zeros((20),dtype=np.int32) 1878 if method in {3,6}: 1879 method_options[0:2] = -1 1880 1881 mi = grid_def_in.npoints 1882 mo = grid_def_out.npoints 1883 1884 # Adjust shape of input array(s) 1885 a, newshp = _adjust_array_shape_for_interp(a, grid_def_in, grid_def_out) 1886 1887 # Call interpolation subroutines according to type of a. 1888 if isinstance(a,np.ndarray): 1889 # Scalar 1890 km = a.shape[0] 1891 if np.any(np.isnan(a)): 1892 ibi = np.ones((km), dtype=np.int32) 1893 li = np.where(np.isnan(a),0,1).astype(np.uint8) 1894 else: 1895 ibi = np.zeros((km), dtype=np.int32) 1896 li = np.zeros(a.shape,dtype=np.uint8) 1897 no, rlat, rlon, ibo, lo, go, iret = iplib.interpolate_scalar( 1898 method, 1899 method_options, 1900 grid_def_in.gdtn, 1901 np.array(grid_def_in.gdt, dtype=np.int32), 1902 grid_def_out.gdtn, 1903 np.array(grid_def_out.gdt, dtype=np.int32), 1904 mi, 1905 mo, 1906 km, 1907 ibi, 1908 li, 1909 a, 1910 ) 1911 out = np.where(lo==0,np.nan,go).reshape(newshp) 1912 elif isinstance(a,tuple): 1913 # Vector 1914 km = a[0].shape[0] 1915 if np.any(np.isnan(a)): 1916 ibi = np.ones((km), dtype=np.int32) 1917 li = np.where(np.isnan(a),0,1).astype(np.uint8) 1918 else: 1919 ibi = np.zeros((km), dtype=np.int32) 1920 li = np.zeros(a[0].shape,dtype=np.uint8) 1921 no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret = iplib.interpolate_vector( 1922 method, 1923 method_options, 1924 grid_def_in.gdtn, 1925 np.array(grid_def_in.gdt, dtype=np.int32), 1926 grid_def_out.gdtn, 1927 np.array(grid_def_out.gdt, dtype=np.int32), 1928 mi, 1929 mo, 1930 km, 1931 ibi, 1932 li, 1933 a[0], 1934 a[1], 1935 ) 1936 uo = np.where(lo==0,np.nan,uo).reshape(newshp) 1937 vo = np.where(lo==0,np.nan,vo).reshape(newshp) 1938 out = (uo,vo) 1939 1940 try: 1941 iplib.openmp_set_num_threads(prev_num_threads) 1942 except(AttributeError): 1943 pass 1944 1945 return out
This is the module-level interpolation function.
This interfaces with the 4-byte (32-bit float) NCEPLIBS-ip library through grib2io's internal iplib Cython extension module.
Parameters
a (numpy.ndarray or tuple): Input data. If
a
is anumpy.ndarray
, scalar interpolation will be performed. Ifa
is atuple
, then vector interpolation will be performed with the assumption that u = a[0] and v = a[1] and are bothnumpy.ndarray
.These data are expected to be in 2-dimensional form with shape (ny, nx) or 3-dimensional (:, ny, nx) where the 1st dimension represents another spatial, temporal, or classification (i.e. ensemble members) dimension. The function will properly flatten the (ny,nx) dimensions into (nx * ny) acceptable for input into the interpolation subroutines. If needed, these data will be converted to
np.float32
.- method: Interpolate method to use. This can either be an integer or string using the following mapping:
Interpolate Scheme | Integer Value |
---|---|
'bilinear' | 0 |
'bicubic' | 1 |
'neighbor' | 2 |
'budget' | 3 |
'spectral' | 4 |
'neighbor-budget' | 6 |
- grid_def_in (grib2io.Grib2GridDef): Grib2GridDef object for the input grid.
- grid_def_out (grib2io.Grib2GridDef): Grib2GridDef object for the output grid or station points.
- method_options (list of ints, optional): Interpolation options. See the NCEPLIBS-ip documentation for more information on how these are used.
- num_threads (int, optional): Number of OpenMP threads to use for interpolation. The default value is 1. If NCEPLIBS-ip and grib2io's iplib extension module was not built with OpenMP, then this keyword argument and value will have no impact.
Returns
- interpolate: Returns a
numpy.ndarray
of dtypenp.float32
when scalar interpolation is performed or atuple
ofnumpy.ndarray
s when vector interpolation is performed with the assumptions that 0-index is the interpolated u and 1-index is the interpolated v.
1948def interpolate_to_stations( 1949 a, 1950 method: Union[int, str], 1951 grid_def_in, 1952 lats, 1953 lons, 1954 method_options=None, 1955 num_threads=1 1956): 1957 """ 1958 Module-level interpolation function for interpolation to stations. 1959 1960 Interfaces with the 4-byte (32-bit float) [NCEPLIBS-ip library](https://github.com/NOAA-EMC/NCEPLIBS-ip) 1961 via grib2io's iplib Cython exntension module. It supports scalar and 1962 vector interpolation according to the type of object a. 1963 1964 Parameters 1965 ---------- 1966 a : numpy.ndarray or tuple 1967 Input data. If `a` is a `numpy.ndarray`, scalar interpolation will be 1968 performed. If `a` is a `tuple`, then vector interpolation will be 1969 performed with the assumption that u = a[0] and v = a[1] and are both 1970 `numpy.ndarray`. 1971 1972 These data are expected to be in 2-dimensional form with shape (ny, nx) 1973 or 3-dimensional (:, ny, nx) where the 1st dimension represents another 1974 spatial, temporal, or classification (i.e. ensemble members) dimension. 1975 The function will properly flatten the (ny,nx) dimensions into (nx * ny) 1976 acceptable for input into the interpolation subroutines. If needed, these 1977 data will be converted to `np.float32`. 1978 method 1979 Interpolate method to use. This can either be an integer or string using 1980 the following mapping: 1981 1982 | Interpolate Scheme | Integer Value | 1983 | :---: | :---: | 1984 | 'bilinear' | 0 | 1985 | 'bicubic' | 1 | 1986 | 'neighbor' | 2 | 1987 | 'budget' | 3 | 1988 | 'spectral' | 4 | 1989 | 'neighbor-budget' | 6 | 1990 1991 grid_def_in : grib2io.Grib2GridDef 1992 Grib2GridDef object for the input grid. 1993 lats : numpy.ndarray or list 1994 Latitudes for station points 1995 lons : numpy.ndarray or list 1996 Longitudes for station points 1997 method_options : list of ints, optional 1998 Interpolation options. See the NCEPLIBS-ip documentation for 1999 more information on how these are used. 2000 num_threads : int, optional 2001 Number of OpenMP threads to use for interpolation. The default 2002 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 2003 was not built with OpenMP, then this keyword argument and value 2004 will have no impact. 2005 2006 Returns 2007 ------- 2008 interpolate_to_stations 2009 Returns a `numpy.ndarray` of dtype `np.float32` when scalar 2010 interpolation is performed or a `tuple` of `numpy.ndarray`s 2011 when vector interpolation is performed with the assumptions 2012 that 0-index is the interpolated u and 1-index is the 2013 interpolated v. 2014 """ 2015 from . import iplib 2016 2017 prev_num_threads = 1 2018 try: 2019 prev_num_threads = iplib.openmp_get_num_threads() 2020 iplib.openmp_set_num_threads(num_threads) 2021 except(AttributeError): 2022 pass 2023 2024 if isinstance(method,int) and method not in _interp_schemes.values(): 2025 raise ValueError('Invalid interpolation method.') 2026 elif isinstance(method,str): 2027 if method in _interp_schemes.keys(): 2028 method = _interp_schemes[method] 2029 else: 2030 raise ValueError('Invalid interpolation method.') 2031 2032 if method_options is None: 2033 method_options = np.zeros((20),dtype=np.int32) 2034 if method in {3,6}: 2035 method_options[0:2] = -1 2036 2037 # Check lats and lons 2038 if isinstance(lats,list): 2039 nlats = len(lats) 2040 elif isinstance(lats,np.ndarray) and len(lats.shape) == 1: 2041 nlats = lats.shape[0] 2042 else: 2043 raise ValueError('Station latitudes must be a list or 1-D NumPy array.') 2044 if isinstance(lons,list): 2045 nlons = len(lons) 2046 elif isinstance(lons,np.ndarray) and len(lons.shape) == 1: 2047 nlons = lons.shape[0] 2048 else: 2049 raise ValueError('Station longitudes must be a list or 1-D NumPy array.') 2050 if nlats != nlons: 2051 raise ValueError('Station lats and lons must be same size.') 2052 2053 mi = grid_def_in.npoints 2054 mo = nlats 2055 2056 # Adjust shape of input array(s) 2057 a, newshp = _adjust_array_shape_for_interp_stations(a, grid_def_in, mo) 2058 2059 # Use gdtn = -1 for stations and an empty template array 2060 gdtn = -1 2061 gdt = np.zeros((200),dtype=np.int32) 2062 2063 # Call interpolation subroutines according to type of a. 2064 if isinstance(a,np.ndarray): 2065 # Scalar 2066 km = a.shape[0] 2067 ibi = np.zeros((km), dtype=np.int32) 2068 li = np.zeros(a.shape,dtype=np.uint8) 2069 no, rlat, rlon, ibo, lo, go, iret = iplib.interpolate_scalar( 2070 method, 2071 method_options, 2072 grid_def_in.gdtn, 2073 np.array(grid_def_in.gdt, dtype=np.int32), 2074 gdtn, 2075 gdt, 2076 mi, 2077 mo, 2078 km, 2079 ibi, 2080 li, 2081 a, 2082 lats=np.array(lats, dtype=np.float32), 2083 lons=np.array(lons, dtype=np.float32), 2084 ) 2085 out = go.reshape(newshp) 2086 2087 elif isinstance(a,tuple): 2088 # Vector 2089 km = a[0].shape[0] 2090 ibi = np.zeros((km), dtype=np.int32) 2091 li = np.zeros(a[0].shape,dtype=np.uint8) 2092 no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret = iplib.interpolate_vector( 2093 method, 2094 method_options, 2095 grid_def_in.gdtn, 2096 np.array(grid_def_in.gdt, dtype=np.int32), 2097 gdtn, 2098 gdt, 2099 mi, 2100 mo, 2101 km, 2102 ibi, 2103 li, 2104 a[0], 2105 a[1], 2106 lats=np.array(lats, dtype=np.float32), 2107 lons=np.array(lons, dtype=np.float32), 2108 ) 2109 out = (uo.reshape(newshp),vo.reshape(newshp)) 2110 2111 try: 2112 iplib.openmp_set_num_threads(prev_num_threads) 2113 except(AttributeError): 2114 pass 2115 2116 return out
Module-level interpolation function for interpolation to stations.
Interfaces with the 4-byte (32-bit float) NCEPLIBS-ip library via grib2io's iplib Cython exntension module. It supports scalar and vector interpolation according to the type of object a.
Parameters
a (numpy.ndarray or tuple): Input data. If
a
is anumpy.ndarray
, scalar interpolation will be performed. Ifa
is atuple
, then vector interpolation will be performed with the assumption that u = a[0] and v = a[1] and are bothnumpy.ndarray
.These data are expected to be in 2-dimensional form with shape (ny, nx) or 3-dimensional (:, ny, nx) where the 1st dimension represents another spatial, temporal, or classification (i.e. ensemble members) dimension. The function will properly flatten the (ny,nx) dimensions into (nx * ny) acceptable for input into the interpolation subroutines. If needed, these data will be converted to
np.float32
.- method: Interpolate method to use. This can either be an integer or string using the following mapping:
Interpolate Scheme | Integer Value |
---|---|
'bilinear' | 0 |
'bicubic' | 1 |
'neighbor' | 2 |
'budget' | 3 |
'spectral' | 4 |
'neighbor-budget' | 6 |
- grid_def_in (grib2io.Grib2GridDef): Grib2GridDef object for the input grid.
- lats (numpy.ndarray or list): Latitudes for station points
- lons (numpy.ndarray or list): Longitudes for station points
- method_options (list of ints, optional): Interpolation options. See the NCEPLIBS-ip documentation for more information on how these are used.
- num_threads (int, optional): Number of OpenMP threads to use for interpolation. The default value is 1. If NCEPLIBS-ip and grib2io's iplib extension module was not built with OpenMP, then this keyword argument and value will have no impact.
Returns
- interpolate_to_stations: Returns a
numpy.ndarray
of dtypenp.float32
when scalar interpolation is performed or atuple
ofnumpy.ndarray
s when vector interpolation is performed with the assumptions that 0-index is the interpolated u and 1-index is the interpolated v.
683class Grib2Message: 684 """ 685 Creation class for a GRIB2 message. 686 687 This class returns a dynamically-created Grib2Message object that 688 inherits from `_Grib2Message` and grid, product, data representation 689 template classes according to the template numbers for the respective 690 sections. If `section3`, `section4`, or `section5` are omitted, then 691 the appropriate keyword arguments for the template number `gdtn=`, 692 `pdtn=`, or `drtn=` must be provided. 693 694 Parameters 695 ---------- 696 section0 697 GRIB2 section 0 array. 698 section1 699 GRIB2 section 1 array. 700 section2 701 Local Use section data. 702 section3 703 GRIB2 section 3 array. 704 section4 705 GRIB2 section 4 array. 706 section5 707 GRIB2 section 5 array. 708 709 Returns 710 ------- 711 Msg 712 A dynamically-create Grib2Message object that inherits from 713 _Grib2Message, a grid definition template class, product 714 definition template class, and a data representation template 715 class. 716 """ 717 def __new__(self, section0: NDArray = np.array([struct.unpack('>I',b'GRIB')[0],0,0,2,0]), 718 section1: NDArray = np.zeros((13),dtype=np.int64), 719 section2: Optional[bytes] = None, 720 section3: Optional[NDArray] = None, 721 section4: Optional[NDArray] = None, 722 section5: Optional[NDArray] = None, *args, **kwargs): 723 724 if np.all(section1==0): 725 try: 726 # Python >= 3.10 727 section1[5:11] = datetime.datetime.fromtimestamp(0, datetime.UTC).timetuple()[:6] 728 except(AttributeError): 729 # Python < 3.10 730 section1[5:11] = datetime.datetime.utcfromtimestamp(0).timetuple()[:6] 731 732 bases = list() 733 if section3 is None: 734 if 'gdtn' in kwargs.keys(): 735 gdtn = kwargs['gdtn'] 736 Gdt = templates.gdt_class_by_gdtn(gdtn) 737 bases.append(Gdt) 738 section3 = np.zeros((Gdt._len+5),dtype=np.int64) 739 section3[4] = gdtn 740 else: 741 raise ValueError("Must provide GRIB2 Grid Definition Template Number or section 3 array") 742 else: 743 gdtn = section3[4] 744 Gdt = templates.gdt_class_by_gdtn(gdtn) 745 bases.append(Gdt) 746 747 if section4 is None: 748 if 'pdtn' in kwargs.keys(): 749 pdtn = kwargs['pdtn'] 750 Pdt = templates.pdt_class_by_pdtn(pdtn) 751 bases.append(Pdt) 752 section4 = np.zeros((Pdt._len+2),dtype=np.int64) 753 section4[1] = pdtn 754 else: 755 raise ValueError("Must provide GRIB2 Production Definition Template Number or section 4 array") 756 else: 757 pdtn = section4[1] 758 Pdt = templates.pdt_class_by_pdtn(pdtn) 759 bases.append(Pdt) 760 761 if section5 is None: 762 if 'drtn' in kwargs.keys(): 763 drtn = kwargs['drtn'] 764 Drt = templates.drt_class_by_drtn(drtn) 765 bases.append(Drt) 766 section5 = np.zeros((Drt._len+2),dtype=np.int64) 767 section5[1] = drtn 768 else: 769 raise ValueError("Must provide GRIB2 Data Representation Template Number or section 5 array") 770 else: 771 drtn = section5[1] 772 Drt = templates.drt_class_by_drtn(drtn) 773 bases.append(Drt) 774 775 # attempt to use existing Msg class if it has already been made with gdtn,pdtn,drtn combo 776 try: 777 Msg = _msg_class_store[f"{gdtn}:{pdtn}:{drtn}"] 778 except KeyError: 779 @dataclass(init=False, repr=False) 780 class Msg(_Grib2Message, *bases): 781 pass 782 _msg_class_store[f"{gdtn}:{pdtn}:{drtn}"] = Msg 783 784 return Msg(section0, section1, section2, section3, section4, section5, *args)
Creation class for a GRIB2 message.
This class returns a dynamically-created Grib2Message object that
inherits from _Grib2Message
and grid, product, data representation
template classes according to the template numbers for the respective
sections. If section3
, section4
, or section5
are omitted, then
the appropriate keyword arguments for the template number gdtn=
,
pdtn=
, or drtn=
must be provided.
Parameters
- section0: GRIB2 section 0 array.
- section1: GRIB2 section 1 array.
- section2: Local Use section data.
- section3: GRIB2 section 3 array.
- section4: GRIB2 section 4 array.
- section5: GRIB2 section 5 array.
Returns
- Msg: A dynamically-create Grib2Message object that inherits from _Grib2Message, a grid definition template class, product definition template class, and a data representation template class.
787@dataclass 788class _Grib2Message: 789 """ 790 GRIB2 Message base class. 791 """ 792 # GRIB2 Sections 793 section0: NDArray = field(init=True,repr=False) 794 section1: NDArray = field(init=True,repr=False) 795 section2: bytes = field(init=True,repr=False) 796 section3: NDArray = field(init=True,repr=False) 797 section4: NDArray = field(init=True,repr=False) 798 section5: NDArray = field(init=True,repr=False) 799 bitMapFlag: templates.Grib2Metadata = field(init=True,repr=False,default=255) 800 801 # Section 0 looked up attributes 802 indicatorSection: NDArray = field(init=False,repr=False,default=templates.IndicatorSection()) 803 discipline: templates.Grib2Metadata = field(init=False,repr=False,default=templates.Discipline()) 804 805 # Section 1 looked up attributes 806 identificationSection: NDArray = field(init=False,repr=False,default=templates.IdentificationSection()) 807 originatingCenter: templates.Grib2Metadata = field(init=False,repr=False,default=templates.OriginatingCenter()) 808 originatingSubCenter: templates.Grib2Metadata = field(init=False,repr=False,default=templates.OriginatingSubCenter()) 809 masterTableInfo: templates.Grib2Metadata = field(init=False,repr=False,default=templates.MasterTableInfo()) 810 localTableInfo: templates.Grib2Metadata = field(init=False,repr=False,default=templates.LocalTableInfo()) 811 significanceOfReferenceTime: templates.Grib2Metadata = field(init=False,repr=False,default=templates.SignificanceOfReferenceTime()) 812 year: int = field(init=False,repr=False,default=templates.Year()) 813 month: int = field(init=False,repr=False,default=templates.Month()) 814 day: int = field(init=False,repr=False,default=templates.Day()) 815 hour: int = field(init=False,repr=False,default=templates.Hour()) 816 minute: int = field(init=False,repr=False,default=templates.Minute()) 817 second: int = field(init=False,repr=False,default=templates.Second()) 818 refDate: datetime.datetime = field(init=False,repr=False,default=templates.RefDate()) 819 productionStatus: templates.Grib2Metadata = field(init=False,repr=False,default=templates.ProductionStatus()) 820 typeOfData: templates.Grib2Metadata = field(init=False,repr=False,default=templates.TypeOfData()) 821 822 # Section 3 looked up common attributes. Other looked up attributes are available according 823 # to the Grid Definition Template. 824 gridDefinitionSection: NDArray = field(init=False,repr=False,default=templates.GridDefinitionSection()) 825 sourceOfGridDefinition: int = field(init=False,repr=False,default=templates.SourceOfGridDefinition()) 826 numberOfDataPoints: int = field(init=False,repr=False,default=templates.NumberOfDataPoints()) 827 interpretationOfListOfNumbers: templates.Grib2Metadata = field(init=False,repr=False,default=templates.InterpretationOfListOfNumbers()) 828 gridDefinitionTemplateNumber: templates.Grib2Metadata = field(init=False,repr=False,default=templates.GridDefinitionTemplateNumber()) 829 gridDefinitionTemplate: list = field(init=False,repr=False,default=templates.GridDefinitionTemplate()) 830 _earthparams: dict = field(init=False,repr=False,default=templates.EarthParams()) 831 _dxsign: float = field(init=False,repr=False,default=templates.DxSign()) 832 _dysign: float = field(init=False,repr=False,default=templates.DySign()) 833 _llscalefactor: float = field(init=False,repr=False,default=templates.LLScaleFactor()) 834 _lldivisor: float = field(init=False,repr=False,default=templates.LLDivisor()) 835 _xydivisor: float = field(init=False,repr=False,default=templates.XYDivisor()) 836 shapeOfEarth: templates.Grib2Metadata = field(init=False,repr=False,default=templates.ShapeOfEarth()) 837 earthShape: str = field(init=False,repr=False,default=templates.EarthShape()) 838 earthRadius: float = field(init=False,repr=False,default=templates.EarthRadius()) 839 earthMajorAxis: float = field(init=False,repr=False,default=templates.EarthMajorAxis()) 840 earthMinorAxis: float = field(init=False,repr=False,default=templates.EarthMinorAxis()) 841 resolutionAndComponentFlags: list = field(init=False,repr=False,default=templates.ResolutionAndComponentFlags()) 842 ny: int = field(init=False,repr=False,default=templates.Ny()) 843 nx: int = field(init=False,repr=False,default=templates.Nx()) 844 scanModeFlags: list = field(init=False,repr=False,default=templates.ScanModeFlags()) 845 projParameters: dict = field(init=False,repr=False,default=templates.ProjParameters()) 846 847 # Section 4 848 productDefinitionTemplateNumber: templates.Grib2Metadata = field(init=False,repr=False,default=templates.ProductDefinitionTemplateNumber()) 849 productDefinitionTemplate: NDArray = field(init=False,repr=False,default=templates.ProductDefinitionTemplate()) 850 851 # Section 5 looked up common attributes. Other looked up attributes are 852 # available according to the Data Representation Template. 853 numberOfPackedValues: int = field(init=False,repr=False,default=templates.NumberOfPackedValues()) 854 dataRepresentationTemplateNumber: templates.Grib2Metadata = field(init=False,repr=False,default=templates.DataRepresentationTemplateNumber()) 855 dataRepresentationTemplate: list = field(init=False,repr=False,default=templates.DataRepresentationTemplate()) 856 typeOfValues: templates.Grib2Metadata = field(init=False,repr=False,default=templates.TypeOfValues()) 857 858 def __post_init__(self): 859 """Set some attributes after init.""" 860 self._auto_nans = _AUTO_NANS 861 self._coordlist = np.zeros((0), dtype=np.float32) 862 self._data = None 863 self._deflist = np.zeros((0), dtype=np.int64) 864 self._msgnum = -1 865 self._ondiskarray = None 866 self._orig_section5 = np.copy(self.section5) 867 self._signature = self._generate_signature() 868 try: 869 self._sha1_section3 = hashlib.sha1(self.section3).hexdigest() 870 except(TypeError): 871 pass 872 self.bitMapFlag = templates.Grib2Metadata(self.bitMapFlag,table='6.0') 873 self.bitmap = None 874 875 @property 876 def _isNDFD(self): 877 """Check if GRIB2 message is from NWS NDFD""" 878 return np.all(self.section1[0:2]==[8,65535]) 879 880 @property 881 def _isAerosol(self): 882 """Check if GRIB2 message contains aerosol data""" 883 is_aero_template = self.productDefinitionTemplateNumber.value in tables.AEROSOL_PDTNS 884 is_aero_param = ((str(self.parameterCategory) == '13') | 885 (str(self.parameterCategory) == '20')) and str(self.parameterNumber) in tables.AEROSOL_PARAMS 886 # Check table 4.205 aerosol presence 887 is_aero_type = (str(self.parameterCategory) == '205' and 888 str(self.parameterNumber) == '1') 889 return is_aero_template or is_aero_param or is_aero_type 890 891 @property 892 def gdtn(self): 893 """Return Grid Definition Template Number""" 894 return self.section3[4] 895 896 897 @property 898 def gdt(self): 899 """Return Grid Definition Template.""" 900 return self.gridDefinitionTemplate 901 902 903 @property 904 def pdtn(self): 905 """Return Product Definition Template Number.""" 906 return self.section4[1] 907 908 909 @property 910 def pdt(self): 911 """Return Product Definition Template.""" 912 return self.productDefinitionTemplate 913 914 915 @property 916 def drtn(self): 917 """Return Data Representation Template Number.""" 918 return self.section5[1] 919 920 921 @property 922 def drt(self): 923 """Return Data Representation Template.""" 924 return self.dataRepresentationTemplate 925 926 927 @property 928 def pdy(self): 929 """Return the PDY ('YYYYMMDD').""" 930 return ''.join([str(i) for i in self.section1[5:8]]) 931 932 933 @property 934 def griddef(self): 935 """Return a Grib2GridDef instance for a GRIB2 message.""" 936 return Grib2GridDef.from_section3(self.section3) 937 938 939 @property 940 def lats(self): 941 """Return grid latitudes.""" 942 return self.latlons()[0] 943 944 945 @property 946 def lons(self): 947 """Return grid longitudes.""" 948 return self.latlons()[1] 949 950 @property 951 def min(self): 952 """Return minimum value of data.""" 953 return np.nanmin(self.data) 954 955 956 @property 957 def max(self): 958 """Return maximum value of data.""" 959 return np.nanmax(self.data) 960 961 962 @property 963 def mean(self): 964 """Return mean value of data.""" 965 return np.nanmean(self.data) 966 967 968 @property 969 def median(self): 970 """Return median value of data.""" 971 return np.nanmedian(self.data) 972 973 974 @property 975 def shape(self): 976 """Return shape of data.""" 977 return self.griddef.shape 978 979 980 def __repr__(self): 981 """ 982 Return an unambiguous string representation of the object. 983 984 Returns 985 ------- 986 repr 987 A string representation of the object, including information from 988 sections 0, 1, 3, 4, 5, and 6. 989 """ 990 info = '' 991 for sect in [0,1,3,4,5,6]: 992 for k,v in self.attrs_by_section(sect,values=True).items(): 993 info += f'Section {sect}: {k} = {v}\n' 994 return info 995 996 def __str__(self): 997 """ 998 Return a readable string representation of the object. 999 1000 Returns 1001 ------- 1002 str 1003 A formatted string representation of the object, including 1004 selected attributes. 1005 """ 1006 return (f'{self._msgnum}:d={self.refDate}:{self.shortName}:' 1007 f'{self.fullName} ({self.units}):{self.level}:' 1008 f'{self.leadTime}') 1009 1010 1011 def _generate_signature(self): 1012 """Generature SHA-1 hash string from GRIB2 integer sections.""" 1013 return hashlib.sha1(np.concatenate((self.section0,self.section1, 1014 self.section3,self.section4, 1015 self.section5))).hexdigest() 1016 1017 1018 def attrs_by_section(self, sect: int, values: bool=False): 1019 """ 1020 Provide a tuple of attribute names for the given GRIB2 section. 1021 1022 Parameters 1023 ---------- 1024 sect 1025 The GRIB2 section number. 1026 values 1027 Optional (default is `False`) argument to return attributes values. 1028 1029 Returns 1030 ------- 1031 attrs_by_section 1032 A list of attribute names or dict of name:value pairs if `values = 1033 True`. 1034 """ 1035 if sect in {0,1,6}: 1036 attrs = templates._section_attrs[sect] 1037 elif sect in {3,4,5}: 1038 def _find_class_index(n): 1039 _key = {3:'Grid', 4:'Product', 5:'Data'} 1040 for i,c in enumerate(self.__class__.__mro__): 1041 if _key[n] in c.__name__: 1042 return i 1043 else: 1044 return [] 1045 if sys.version_info.minor <= 8: 1046 attrs = templates._section_attrs[sect]+\ 1047 [a for a in dir(self.__class__.__mro__[_find_class_index(sect)]) if not a.startswith('_')] 1048 else: 1049 attrs = templates._section_attrs[sect]+\ 1050 self.__class__.__mro__[_find_class_index(sect)]._attrs() 1051 else: 1052 attrs = [] 1053 if values: 1054 return {k:getattr(self,k) for k in attrs} 1055 else: 1056 return attrs 1057 1058 1059 def pack(self): 1060 """ 1061 Pack GRIB2 section data into a binary message. 1062 1063 It is the user's responsibility to populate the GRIB2 section 1064 information with appropriate metadata. 1065 """ 1066 # Create beginning of packed binary message with section 0 and 1 data. 1067 self._sections = [] 1068 self._msg,self._pos = g2clib.grib2_create(self.indicatorSection[2:4],self.identificationSection) 1069 self._sections += [0,1] 1070 1071 # Add section 2 if present. 1072 if isinstance(self.section2,bytes) and len(self.section2) > 0: 1073 self._msg,self._pos = g2clib.grib2_addlocal(self._msg,self.section2) 1074 self._sections.append(2) 1075 1076 # Add section 3. 1077 self.section3[1] = self.nx * self.ny 1078 self._msg,self._pos = g2clib.grib2_addgrid(self._msg,self.gridDefinitionSection, 1079 self.gridDefinitionTemplate, 1080 self._deflist) 1081 self._sections.append(3) 1082 1083 # Prepare data. 1084 if self._data is None: 1085 if self._ondiskarray is None: 1086 raise ValueError("Grib2Message object has no data, thus it cannot be packed.") 1087 field = np.copy(self.data) 1088 if self.scanModeFlags is not None: 1089 if self.scanModeFlags[3]: 1090 fieldsave = field.astype('f') # Casting makes a copy 1091 field[1::2,:] = fieldsave[1::2,::-1] 1092 fld = field.astype('f') 1093 1094 # Prepare bitmap, if necessary 1095 bitmapflag = self.bitMapFlag.value 1096 if bitmapflag == 0: 1097 if self.bitmap is not None: 1098 bmap = np.ravel(self.bitmap).astype(DEFAULT_NUMPY_INT) 1099 else: 1100 bmap = np.ravel(np.where(np.isnan(fld),0,1)).astype(DEFAULT_NUMPY_INT) 1101 else: 1102 bmap = None 1103 1104 # Prepare data for packing if nans are present 1105 fld = np.ravel(fld) 1106 if bitmapflag in {0,254}: 1107 fld = np.where(np.isnan(fld),0,fld) 1108 else: 1109 if np.isnan(fld).any(): 1110 if hasattr(self,'priMissingValue'): 1111 fld = np.where(np.isnan(fld),self.priMissingValue,fld) 1112 if hasattr(self,'_missvalmap'): 1113 if hasattr(self,'priMissingValue'): 1114 fld = np.where(self._missvalmap==1,self.priMissingValue,fld) 1115 if hasattr(self,'secMissingValue'): 1116 fld = np.where(self._missvalmap==2,self.secMissingValue,fld) 1117 1118 # Add sections 4, 5, 6, and 7. 1119 self._msg,self._pos = g2clib.grib2_addfield(self._msg,self.pdtn, 1120 self.productDefinitionTemplate, 1121 self._coordlist, 1122 self.drtn, 1123 self.dataRepresentationTemplate, 1124 fld, 1125 bitmapflag, 1126 bmap) 1127 self._sections.append(4) 1128 self._sections.append(5) 1129 self._sections.append(6) 1130 self._sections.append(7) 1131 1132 # Finalize GRIB2 message with section 8. 1133 self._msg, self._pos = g2clib.grib2_end(self._msg) 1134 self._sections.append(8) 1135 self.section0[-1] = len(self._msg) 1136 1137 1138 @property 1139 def data(self) -> np.array: 1140 """Access the unpacked data values.""" 1141 if self._data is None: 1142 if self._auto_nans != _AUTO_NANS: 1143 self._data = self._ondiskarray 1144 self._data = np.asarray(self._ondiskarray) 1145 return self._data 1146 1147 1148 @data.setter 1149 def data(self, data): 1150 if not isinstance(data, np.ndarray): 1151 raise ValueError('Grib2Message data only supports numpy arrays') 1152 self._data = data 1153 1154 1155 def flush_data(self): 1156 """ 1157 Flush the unpacked data values from the Grib2Message object. 1158 1159 Note: If the Grib2Message object was constructed from "scratch" (i.e. 1160 not read from file), this method will remove the data array from 1161 the object and it cannot be recovered. 1162 """ 1163 self._data = None 1164 self.bitmap = None 1165 1166 1167 def __getitem__(self, item): 1168 return self.data[item] 1169 1170 1171 def __setitem__(self, item): 1172 raise NotImplementedError('assignment of data not supported via setitem') 1173 1174 1175 def latlons(self, *args, **kwrgs): 1176 """Alias for `grib2io.Grib2Message.grid` method.""" 1177 return self.grid(*args, **kwrgs) 1178 1179 1180 def grid(self, unrotate: bool=True): 1181 """ 1182 Return lats,lons (in degrees) of grid. 1183 1184 Currently can handle reg. lat/lon,cglobal Gaussian, mercator, 1185 stereographic, lambert conformal, albers equal-area, space-view and 1186 azimuthal equidistant grids. 1187 1188 Parameters 1189 ---------- 1190 unrotate 1191 If `True` [DEFAULT], and grid is rotated lat/lon, then unrotate the 1192 grid, otherwise `False`, do not. 1193 1194 Returns 1195 ------- 1196 lats, lons : numpy.ndarray 1197 Returns two numpy.ndarrays with dtype=numpy.float32 of grid 1198 latitudes and longitudes in units of degrees. 1199 """ 1200 if self._sha1_section3 in _latlon_datastore.keys(): 1201 return (_latlon_datastore[self._sha1_section3]['latitude'], 1202 _latlon_datastore[self._sha1_section3]['longitude']) 1203 gdtn = self.gridDefinitionTemplateNumber.value 1204 gdtmpl = self.gridDefinitionTemplate 1205 reggrid = self.gridDefinitionSection[2] == 0 # This means regular 2-d grid 1206 if gdtn == 0: 1207 # Regular lat/lon grid 1208 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1209 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1210 dlon = self.gridlengthXDirection 1211 dlat = self.gridlengthYDirection 1212 if lon2 < lon1 and dlon < 0: lon1 = -lon1 1213 lats = np.linspace(lat1,lat2,self.ny) 1214 if reggrid: 1215 lons = np.linspace(lon1,lon2,self.nx) 1216 else: 1217 lons = np.linspace(lon1,lon2,self.ny*2) 1218 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1219 elif gdtn == 1: # Rotated Lat/Lon grid 1220 pj = pyproj.Proj(self.projParameters) 1221 lat1,lon1 = self.latitudeFirstGridpoint,self.longitudeFirstGridpoint 1222 lat2,lon2 = self.latitudeLastGridpoint,self.longitudeLastGridpoint 1223 if lon1 > 180.0: lon1 -= 360.0 1224 if lon2 > 180.0: lon2 -= 360.0 1225 lats = np.linspace(lat1,lat2,self.ny) 1226 lons = np.linspace(lon1,lon2,self.nx) 1227 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1228 if unrotate: 1229 from grib2io.utils import rotated_grid 1230 lats,lons = rotated_grid.unrotate(lats,lons,self.anglePoleRotation, 1231 self.latitudeSouthernPole, 1232 self.longitudeSouthernPole) 1233 elif gdtn == 40: # Gaussian grid (only works for global!) 1234 from grib2io.utils.gauss_grid import gaussian_latitudes 1235 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1236 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1237 nlats = self.ny 1238 if not reggrid: # Reduced Gaussian grid. 1239 nlons = 2*nlats 1240 dlon = 360./nlons 1241 else: 1242 nlons = self.nx 1243 dlon = self.gridlengthXDirection 1244 lons = np.linspace(lon1,lon2,nlons) 1245 # Compute Gaussian lats (north to south) 1246 lats = gaussian_latitudes(nlats) 1247 if lat1 < lat2: # reverse them if necessary 1248 lats = lats[::-1] 1249 lons,lats = np.meshgrid(lons,lats) 1250 elif gdtn in {10,20,30,31,110}: 1251 # Mercator, Lambert Conformal, Stereographic, Albers Equal Area, 1252 # Azimuthal Equidistant 1253 dx,dy = self.gridlengthXDirection, self.gridlengthYDirection 1254 lon1,lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1255 pj = pyproj.Proj(self.projParameters) 1256 llcrnrx, llcrnry = pj(lon1,lat1) 1257 x = llcrnrx+dx*np.arange(self.nx) 1258 y = llcrnry+dy*np.arange(self.ny) 1259 x,y = np.meshgrid(x, y) 1260 lons,lats = pj(x, y, inverse=True) 1261 elif gdtn == 90: 1262 # Satellite Projection 1263 dx = self.gridlengthXDirection 1264 dy = self.gridlengthYDirection 1265 pj = pyproj.Proj(self.projParameters) 1266 x = dx*np.indices((self.ny,self.nx),'f')[1,:,:] 1267 x -= 0.5*x.max() 1268 y = dy*np.indices((self.ny,self.nx),'f')[0,:,:] 1269 y -= 0.5*y.max() 1270 lons,lats = pj(x,y,inverse=True) 1271 # Set lons,lats to 1.e30 where undefined 1272 abslons = np.fabs(lons) 1273 abslats = np.fabs(lats) 1274 lons = np.where(abslons < 1.e20, lons, 1.e30) 1275 lats = np.where(abslats < 1.e20, lats, 1.e30) 1276 elif gdtn == 32769: 1277 # Special NCEP Grid, Rotated Lat/Lon, Arakawa E-Grid (Non-Staggered) 1278 from grib2io.utils import arakawa_rotated_grid 1279 from grib2io.utils.rotated_grid import DEG2RAD 1280 di, dj = 0.0, 0.0 1281 do_180 = False 1282 idir = 1 if self.scanModeFlags[0] == 0 else -1 1283 jdir = -1 if self.scanModeFlags[1] == 0 else 1 1284 grid_oriented = 0 if self.resolutionAndComponentFlags[4] == 0 else 1 1285 do_rot = 1 if grid_oriented == 1 else 0 1286 la1 = self.latitudeFirstGridpoint 1287 lo1 = self.longitudeFirstGridpoint 1288 clon = self.longitudeCenterGridpoint 1289 clat = self.latitudeCenterGridpoint 1290 lasp = clat - 90.0 1291 losp = clon 1292 llat, llon = arakawa_rotated_grid.ll2rot(la1,lo1,lasp,losp) 1293 la2, lo2 = arakawa_rotated_grid.rot2ll(-llat,-llon,lasp,losp) 1294 rlat = -llat 1295 rlon = -llon 1296 if self.nx == 1: 1297 di = 0.0 1298 elif idir == 1: 1299 ti = rlon 1300 while ti < llon: 1301 ti += 360.0 1302 di = (ti - llon)/float(self.nx-1) 1303 else: 1304 ti = llon 1305 while ti < rlon: 1306 ti += 360.0 1307 di = (ti - rlon)/float(self.nx-1) 1308 if self.ny == 1: 1309 dj = 0.0 1310 else: 1311 dj = (rlat - llat)/float(self.ny-1) 1312 if dj < 0.0: 1313 dj = -dj 1314 if idir == 1: 1315 if llon > rlon: 1316 llon -= 360.0 1317 if llon < 0 and rlon > 0: 1318 do_180 = True 1319 else: 1320 if rlon > llon: 1321 rlon -= 360.0 1322 if rlon < 0 and llon > 0: 1323 do_180 = True 1324 xlat1d = llat + (np.arange(self.ny)*jdir*dj) 1325 xlon1d = llon + (np.arange(self.nx)*idir*di) 1326 xlons, xlats = np.meshgrid(xlon1d,xlat1d) 1327 rot2ll_vectorized = np.vectorize(arakawa_rotated_grid.rot2ll) 1328 lats, lons = rot2ll_vectorized(xlats,xlons,lasp,losp) 1329 if do_180: 1330 lons = np.where(lons>180.0,lons-360.0,lons) 1331 vector_rotation_angles_vectorized = np.vectorize(arakawa_rotated_grid.vector_rotation_angles) 1332 rots = vector_rotation_angles_vectorized(lats, lons, clat, losp, xlats) 1333 del xlat1d, xlon1d, xlats, xlons 1334 else: 1335 raise ValueError('Unsupported grid') 1336 1337 _latlon_datastore[self._sha1_section3] = dict(latitude=lats,longitude=lons) 1338 try: 1339 _latlon_datastore[self._sha1_section3]['vector_rotation_angles'] = rots 1340 except(NameError): 1341 pass 1342 1343 return lats, lons 1344 1345 1346 def map_keys(self): 1347 """ 1348 Unpack data grid replacing integer values with strings. 1349 1350 These types of fields are categorical or classifications where data 1351 values do not represent an observable or predictable physical quantity. 1352 An example of such a field would be [Dominant Precipitation Type - 1353 DPTYPE](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-201.shtml) 1354 1355 Returns 1356 ------- 1357 map_keys 1358 numpy.ndarray of string values per element. 1359 """ 1360 hold_auto_nans = _AUTO_NANS 1361 set_auto_nans(False) 1362 1363 if (np.all(self.section1[0:2] == [7, 14]) and self.shortName == "PWTHER") or ( 1364 self._isNDFD and self.shortName in {"WX", "WWA"} 1365 ): 1366 keys = utils.decode_wx_strings(self.section2) 1367 if hasattr(self,'priMissingValue') and self.priMissingValue not in [None,0]: 1368 keys[int(self.priMissingValue)] = 'Missing' 1369 if hasattr(self,'secMissingValue') and self.secMissingValue not in [None,0]: 1370 keys[int(self.secMissingValue)] = 'Missing' 1371 u,inv = np.unique(self.data,return_inverse=True) 1372 fld = np.array([keys[x] for x in u])[inv].reshape(self.data.shape) 1373 else: 1374 # For data whose units are defined in a code table (i.e. classification or mask) 1375 tblname = re.findall(r'\d\.\d+',self.units,re.IGNORECASE)[0] 1376 fld = self.data.astype(np.int32).astype(str) 1377 tbl = tables.get_table(tblname,expand=True) 1378 for val in np.unique(fld): 1379 fld = np.where(fld==val,tbl[val],fld) 1380 set_auto_nans(hold_auto_nans) 1381 return fld 1382 1383 1384 def to_bytes(self, validate: bool=True): 1385 """ 1386 Return packed GRIB2 message in bytes format. 1387 1388 This will be useful for exporting data in non-file formats. For example, 1389 can be used to output grib data directly to S3 using the boto3 client 1390 without the need to write a temporary file to upload first. 1391 1392 Parameters 1393 ---------- 1394 validate: default=True 1395 If `True` (DEFAULT), validates first/last four bytes for proper 1396 formatting, else returns None. If `False`, message is output as is. 1397 1398 Returns 1399 ------- 1400 to_bytes 1401 Returns GRIB2 formatted message as bytes. 1402 """ 1403 if hasattr(self,'_msg'): 1404 if validate: 1405 if self.validate(): 1406 return self._msg 1407 else: 1408 return None 1409 else: 1410 return self._msg 1411 else: 1412 return None 1413 1414 1415 def interpolate(self, method, grid_def_out, method_options=None, drtn=None, 1416 num_threads=1): 1417 """ 1418 Grib2Message Interpolator 1419 1420 Performs spatial interpolation via the [NCEPLIBS-ip 1421 library](https://github.com/NOAA-EMC/NCEPLIBS-ip). This interpolate 1422 method only supports scalar interpolation. If you need to perform 1423 vector interpolation, use the module-level `grib2io.interpolate` 1424 function. 1425 1426 Parameters 1427 ---------- 1428 method 1429 Interpolate method to use. This can either be an integer or string 1430 using the following mapping: 1431 1432 | Interpolate Scheme | Integer Value | 1433 | :---: | :---: | 1434 | 'bilinear' | 0 | 1435 | 'bicubic' | 1 | 1436 | 'neighbor' | 2 | 1437 | 'budget' | 3 | 1438 | 'spectral' | 4 | 1439 | 'neighbor-budget' | 6 | 1440 1441 grid_def_out : grib2io.Grib2GridDef 1442 Grib2GridDef object of the output grid. 1443 method_options : list of ints, optional 1444 Interpolation options. See the NCEPLIBS-ip documentation for 1445 more information on how these are used. 1446 drtn 1447 Data Representation Template to be used for the returned 1448 interpolated GRIB2 message. When `None`, the data representation 1449 template of the source GRIB2 message is used. Once again, it is the 1450 user's responsibility to properly set the Data Representation 1451 Template attributes. 1452 num_threads : int, optional 1453 Number of OpenMP threads to use for interpolation. The default 1454 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 1455 was not built with OpenMP, then this keyword argument and value 1456 will have no impact. 1457 1458 Returns 1459 ------- 1460 interpolate 1461 If interpolating to a grid, a new Grib2Message object is returned. 1462 The GRIB2 metadata of the new Grib2Message object is identical to 1463 the input except where required to be different because of the new 1464 grid specs and possibly a new data representation template. 1465 1466 If interpolating to station points, the interpolated data values are 1467 returned as a numpy.ndarray. 1468 """ 1469 section0 = self.section0 1470 section0[-1] = 0 1471 gds = [0, grid_def_out.npoints, 0, 255, grid_def_out.gdtn] 1472 section3 = np.concatenate((gds,grid_def_out.gdt)) 1473 drtn = self.drtn if drtn is None else drtn 1474 1475 msg = Grib2Message(section0,self.section1,self.section2,section3, 1476 self.section4,None,self.bitMapFlag.value,drtn=drtn) 1477 1478 msg._msgnum = -1 1479 msg._deflist = self._deflist 1480 msg._coordlist = self._coordlist 1481 shape = (msg.ny,msg.nx) 1482 ndim = 2 1483 if msg.typeOfValues == 0: 1484 dtype = 'float32' 1485 elif msg.typeOfValues == 1: 1486 dtype = 'int32' 1487 msg._data = interpolate(self.data,method,Grib2GridDef.from_section3(self.section3),grid_def_out, 1488 method_options=method_options,num_threads=num_threads).reshape(msg.ny,msg.nx) 1489 msg.section5[0] = grid_def_out.npoints 1490 return msg 1491 1492 1493 def subset(self, lats, lons): 1494 """ 1495 Return a spatial subset. 1496 1497 Currently only supports regular grids of the following types: 1498 1499 | Grid Type | gdtn | 1500 | :---: | :---: | 1501 | Latitude/Longitude, Equidistant Cylindrical, or Plate Carree | 0 | 1502 | Rotated Latitude/Longitude | 1 | 1503 | Mercator | 10 | 1504 | Polar Stereographic | 20 | 1505 | Lambert Conformal | 30 | 1506 | Albers Equal-Area | 31 | 1507 | Gaussian Latitude/Longitude | 40 | 1508 | Equatorial Azimuthal Equidistant Projection | 110 | 1509 1510 Parameters 1511 ---------- 1512 lats 1513 List or tuple of latitudes. The minimum and maximum latitudes will 1514 be used to define the southern and northern boundaries. 1515 1516 The order of the latitudes is not important. The function will 1517 determine which is the minimum and maximum. 1518 1519 The latitudes should be in decimal degrees with 0.0 at the equator, 1520 positive values in the northern hemisphere increasing to 90, and 1521 negative values in the southern hemisphere decreasing to -90. 1522 lons 1523 List or tuple of longitudes. The minimum and maximum longitudes 1524 will be used to define the western and eastern boundaries. 1525 1526 The order of the longitudes is not important. The function will 1527 determine which is the minimum and maximum. 1528 1529 GRIB2 longitudes should be in decimal degrees with 0.0 at the prime 1530 meridian, positive values increasing eastward to 360. There are no 1531 negative GRIB2 longitudes. 1532 1533 The typical west longitudes that start at 0.0 at the prime meridian 1534 and decrease to -180 westward, are converted to GRIB2 longitudes by 1535 '360 - (absolute value of the west longitude)' where typical 1536 eastern longitudes are unchanged as GRIB2 longitudes. 1537 1538 Returns 1539 ------- 1540 subset 1541 A spatial subset of a GRIB2 message. 1542 """ 1543 if self.gdtn not in [0, 1, 10, 20, 30, 31, 40, 110]: 1544 raise ValueError( 1545 """ 1546 1547Subset only works for 1548 Latitude/Longitude, Equidistant Cylindrical, or Plate Carree (gdtn=0) 1549 Rotated Latitude/Longitude (gdtn=1) 1550 Mercator (gdtn=10) 1551 Polar Stereographic (gdtn=20) 1552 Lambert Conformal (gdtn=30) 1553 Albers Equal-Area (gdtn=31) 1554 Gaussian Latitude/Longitude (gdtn=40) 1555 Equatorial Azimuthal Equidistant Projection (gdtn=110) 1556 1557""" 1558 ) 1559 1560 if self.nx == 0 or self.ny == 0: 1561 raise ValueError( 1562 """ 1563 1564Subset only works for regular grids. 1565 1566""" 1567 ) 1568 1569 newmsg = Grib2Message( 1570 np.copy(self.section0), 1571 np.copy(self.section1), 1572 np.copy(self.section2), 1573 np.copy(self.section3), 1574 np.copy(self.section4), 1575 np.copy(self.section5), 1576 ) 1577 1578 msglats, msglons = self.grid() 1579 1580 la1 = np.max(lats) 1581 lo1 = np.min(lons) 1582 la2 = np.min(lats) 1583 lo2 = np.max(lons) 1584 1585 # Find the indices of the first and last grid points to the nearest 1586 # lat/lon values in the grid. 1587 first_lat = np.abs(msglats - la1) 1588 first_lon = np.abs(msglons - lo1) 1589 max_idx = np.maximum(first_lat, first_lon) 1590 first_j, first_i = np.where(max_idx == np.min(max_idx)) 1591 1592 last_lat = np.abs(msglats - la2) 1593 last_lon = np.abs(msglons - lo2) 1594 max_idx = np.maximum(last_lat, last_lon) 1595 last_j, last_i = np.where(max_idx == np.min(max_idx)) 1596 1597 setattr(newmsg, "latitudeFirstGridpoint", msglats[first_j[0], first_i[0]]) 1598 setattr(newmsg, "longitudeFirstGridpoint", msglons[first_j[0], first_i[0]]) 1599 setattr(newmsg, "nx", np.abs(first_i[0] - last_i[0])) 1600 setattr(newmsg, "ny", np.abs(first_j[0] - last_j[0])) 1601 1602 # Set *LastGridpoint attributes even if only used for gdtn=[0, 1, 40]. 1603 # This information is used to subset xarray datasets and even though 1604 # unnecessary for some supported grid types, it won't affect a grib2io 1605 # message to set them. 1606 setattr(newmsg, "latitudeLastGridpoint", msglats[last_j[0], last_i[0]]) 1607 setattr(newmsg, "longitudeLastGridpoint", msglons[last_j[0], last_i[0]]) 1608 1609 setattr( 1610 newmsg, 1611 "data", 1612 self.data[ 1613 min(first_j[0], last_j[0]) : max(first_j[0], last_j[0]), 1614 min(first_i[0], last_i[0]) : max(first_i[0], last_i[0]), 1615 ].copy(), 1616 ) 1617 1618 # Need to set the newmsg._sha1_section3 to a blank string so the grid 1619 # method ignores the cached lat/lon values. This will force the grid 1620 # method to recompute the lat/lon values for the subsetted grid. 1621 newmsg._sha1_section3 = "" 1622 newmsg.grid() 1623 1624 return newmsg 1625 1626 def validate(self): 1627 """ 1628 Validate a complete GRIB2 message. 1629 1630 The g2c library does its own internal validation when g2_gribend() is called, but 1631 we will check in grib2io also. The validation checks if the first 4 bytes in 1632 self._msg is 'GRIB' and '7777' as the last 4 bytes and that the message length in 1633 section 0 equals the length of the packed message. 1634 1635 Returns 1636 ------- 1637 `True` if the packed GRIB2 message is complete and well-formed, `False` otherwise. 1638 """ 1639 valid = False 1640 if hasattr(self,'_msg'): 1641 if self._msg[0:4]+self._msg[-4:] == b'GRIB7777': 1642 if self.section0[-1] == len(self._msg): 1643 valid = True 1644 return valid
GRIB2 Message base class.
GRIB2 Section 1, Identification Section
GRIB2 Section 3, Grid Definition Section
Product Definition Template
891 @property 892 def gdtn(self): 893 """Return Grid Definition Template Number""" 894 return self.section3[4]
Return Grid Definition Template Number
897 @property 898 def gdt(self): 899 """Return Grid Definition Template.""" 900 return self.gridDefinitionTemplate
Return Grid Definition Template.
903 @property 904 def pdtn(self): 905 """Return Product Definition Template Number.""" 906 return self.section4[1]
Return Product Definition Template Number.
909 @property 910 def pdt(self): 911 """Return Product Definition Template.""" 912 return self.productDefinitionTemplate
Return Product Definition Template.
915 @property 916 def drtn(self): 917 """Return Data Representation Template Number.""" 918 return self.section5[1]
Return Data Representation Template Number.
921 @property 922 def drt(self): 923 """Return Data Representation Template.""" 924 return self.dataRepresentationTemplate
Return Data Representation Template.
927 @property 928 def pdy(self): 929 """Return the PDY ('YYYYMMDD').""" 930 return ''.join([str(i) for i in self.section1[5:8]])
Return the PDY ('YYYYMMDD').
933 @property 934 def griddef(self): 935 """Return a Grib2GridDef instance for a GRIB2 message.""" 936 return Grib2GridDef.from_section3(self.section3)
Return a Grib2GridDef instance for a GRIB2 message.
950 @property 951 def min(self): 952 """Return minimum value of data.""" 953 return np.nanmin(self.data)
Return minimum value of data.
956 @property 957 def max(self): 958 """Return maximum value of data.""" 959 return np.nanmax(self.data)
Return maximum value of data.
962 @property 963 def mean(self): 964 """Return mean value of data.""" 965 return np.nanmean(self.data)
Return mean value of data.
968 @property 969 def median(self): 970 """Return median value of data.""" 971 return np.nanmedian(self.data)
Return median value of data.
1018 def attrs_by_section(self, sect: int, values: bool=False): 1019 """ 1020 Provide a tuple of attribute names for the given GRIB2 section. 1021 1022 Parameters 1023 ---------- 1024 sect 1025 The GRIB2 section number. 1026 values 1027 Optional (default is `False`) argument to return attributes values. 1028 1029 Returns 1030 ------- 1031 attrs_by_section 1032 A list of attribute names or dict of name:value pairs if `values = 1033 True`. 1034 """ 1035 if sect in {0,1,6}: 1036 attrs = templates._section_attrs[sect] 1037 elif sect in {3,4,5}: 1038 def _find_class_index(n): 1039 _key = {3:'Grid', 4:'Product', 5:'Data'} 1040 for i,c in enumerate(self.__class__.__mro__): 1041 if _key[n] in c.__name__: 1042 return i 1043 else: 1044 return [] 1045 if sys.version_info.minor <= 8: 1046 attrs = templates._section_attrs[sect]+\ 1047 [a for a in dir(self.__class__.__mro__[_find_class_index(sect)]) if not a.startswith('_')] 1048 else: 1049 attrs = templates._section_attrs[sect]+\ 1050 self.__class__.__mro__[_find_class_index(sect)]._attrs() 1051 else: 1052 attrs = [] 1053 if values: 1054 return {k:getattr(self,k) for k in attrs} 1055 else: 1056 return attrs
Provide a tuple of attribute names for the given GRIB2 section.
Parameters
- sect: The GRIB2 section number.
- values: Optional (default is
False
) argument to return attributes values.
Returns
- attrs_by_section: A list of attribute names or dict of name:value pairs if
values = True
.
1059 def pack(self): 1060 """ 1061 Pack GRIB2 section data into a binary message. 1062 1063 It is the user's responsibility to populate the GRIB2 section 1064 information with appropriate metadata. 1065 """ 1066 # Create beginning of packed binary message with section 0 and 1 data. 1067 self._sections = [] 1068 self._msg,self._pos = g2clib.grib2_create(self.indicatorSection[2:4],self.identificationSection) 1069 self._sections += [0,1] 1070 1071 # Add section 2 if present. 1072 if isinstance(self.section2,bytes) and len(self.section2) > 0: 1073 self._msg,self._pos = g2clib.grib2_addlocal(self._msg,self.section2) 1074 self._sections.append(2) 1075 1076 # Add section 3. 1077 self.section3[1] = self.nx * self.ny 1078 self._msg,self._pos = g2clib.grib2_addgrid(self._msg,self.gridDefinitionSection, 1079 self.gridDefinitionTemplate, 1080 self._deflist) 1081 self._sections.append(3) 1082 1083 # Prepare data. 1084 if self._data is None: 1085 if self._ondiskarray is None: 1086 raise ValueError("Grib2Message object has no data, thus it cannot be packed.") 1087 field = np.copy(self.data) 1088 if self.scanModeFlags is not None: 1089 if self.scanModeFlags[3]: 1090 fieldsave = field.astype('f') # Casting makes a copy 1091 field[1::2,:] = fieldsave[1::2,::-1] 1092 fld = field.astype('f') 1093 1094 # Prepare bitmap, if necessary 1095 bitmapflag = self.bitMapFlag.value 1096 if bitmapflag == 0: 1097 if self.bitmap is not None: 1098 bmap = np.ravel(self.bitmap).astype(DEFAULT_NUMPY_INT) 1099 else: 1100 bmap = np.ravel(np.where(np.isnan(fld),0,1)).astype(DEFAULT_NUMPY_INT) 1101 else: 1102 bmap = None 1103 1104 # Prepare data for packing if nans are present 1105 fld = np.ravel(fld) 1106 if bitmapflag in {0,254}: 1107 fld = np.where(np.isnan(fld),0,fld) 1108 else: 1109 if np.isnan(fld).any(): 1110 if hasattr(self,'priMissingValue'): 1111 fld = np.where(np.isnan(fld),self.priMissingValue,fld) 1112 if hasattr(self,'_missvalmap'): 1113 if hasattr(self,'priMissingValue'): 1114 fld = np.where(self._missvalmap==1,self.priMissingValue,fld) 1115 if hasattr(self,'secMissingValue'): 1116 fld = np.where(self._missvalmap==2,self.secMissingValue,fld) 1117 1118 # Add sections 4, 5, 6, and 7. 1119 self._msg,self._pos = g2clib.grib2_addfield(self._msg,self.pdtn, 1120 self.productDefinitionTemplate, 1121 self._coordlist, 1122 self.drtn, 1123 self.dataRepresentationTemplate, 1124 fld, 1125 bitmapflag, 1126 bmap) 1127 self._sections.append(4) 1128 self._sections.append(5) 1129 self._sections.append(6) 1130 self._sections.append(7) 1131 1132 # Finalize GRIB2 message with section 8. 1133 self._msg, self._pos = g2clib.grib2_end(self._msg) 1134 self._sections.append(8) 1135 self.section0[-1] = len(self._msg)
Pack GRIB2 section data into a binary message.
It is the user's responsibility to populate the GRIB2 section information with appropriate metadata.
1138 @property 1139 def data(self) -> np.array: 1140 """Access the unpacked data values.""" 1141 if self._data is None: 1142 if self._auto_nans != _AUTO_NANS: 1143 self._data = self._ondiskarray 1144 self._data = np.asarray(self._ondiskarray) 1145 return self._data
Access the unpacked data values.
1155 def flush_data(self): 1156 """ 1157 Flush the unpacked data values from the Grib2Message object. 1158 1159 Note: If the Grib2Message object was constructed from "scratch" (i.e. 1160 not read from file), this method will remove the data array from 1161 the object and it cannot be recovered. 1162 """ 1163 self._data = None 1164 self.bitmap = None
Flush the unpacked data values from the Grib2Message object.
Note: If the Grib2Message object was constructed from "scratch" (i.e. not read from file), this method will remove the data array from the object and it cannot be recovered.
1175 def latlons(self, *args, **kwrgs): 1176 """Alias for `grib2io.Grib2Message.grid` method.""" 1177 return self.grid(*args, **kwrgs)
Alias for grib2io.Grib2Message.grid
method.
1180 def grid(self, unrotate: bool=True): 1181 """ 1182 Return lats,lons (in degrees) of grid. 1183 1184 Currently can handle reg. lat/lon,cglobal Gaussian, mercator, 1185 stereographic, lambert conformal, albers equal-area, space-view and 1186 azimuthal equidistant grids. 1187 1188 Parameters 1189 ---------- 1190 unrotate 1191 If `True` [DEFAULT], and grid is rotated lat/lon, then unrotate the 1192 grid, otherwise `False`, do not. 1193 1194 Returns 1195 ------- 1196 lats, lons : numpy.ndarray 1197 Returns two numpy.ndarrays with dtype=numpy.float32 of grid 1198 latitudes and longitudes in units of degrees. 1199 """ 1200 if self._sha1_section3 in _latlon_datastore.keys(): 1201 return (_latlon_datastore[self._sha1_section3]['latitude'], 1202 _latlon_datastore[self._sha1_section3]['longitude']) 1203 gdtn = self.gridDefinitionTemplateNumber.value 1204 gdtmpl = self.gridDefinitionTemplate 1205 reggrid = self.gridDefinitionSection[2] == 0 # This means regular 2-d grid 1206 if gdtn == 0: 1207 # Regular lat/lon grid 1208 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1209 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1210 dlon = self.gridlengthXDirection 1211 dlat = self.gridlengthYDirection 1212 if lon2 < lon1 and dlon < 0: lon1 = -lon1 1213 lats = np.linspace(lat1,lat2,self.ny) 1214 if reggrid: 1215 lons = np.linspace(lon1,lon2,self.nx) 1216 else: 1217 lons = np.linspace(lon1,lon2,self.ny*2) 1218 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1219 elif gdtn == 1: # Rotated Lat/Lon grid 1220 pj = pyproj.Proj(self.projParameters) 1221 lat1,lon1 = self.latitudeFirstGridpoint,self.longitudeFirstGridpoint 1222 lat2,lon2 = self.latitudeLastGridpoint,self.longitudeLastGridpoint 1223 if lon1 > 180.0: lon1 -= 360.0 1224 if lon2 > 180.0: lon2 -= 360.0 1225 lats = np.linspace(lat1,lat2,self.ny) 1226 lons = np.linspace(lon1,lon2,self.nx) 1227 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1228 if unrotate: 1229 from grib2io.utils import rotated_grid 1230 lats,lons = rotated_grid.unrotate(lats,lons,self.anglePoleRotation, 1231 self.latitudeSouthernPole, 1232 self.longitudeSouthernPole) 1233 elif gdtn == 40: # Gaussian grid (only works for global!) 1234 from grib2io.utils.gauss_grid import gaussian_latitudes 1235 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1236 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1237 nlats = self.ny 1238 if not reggrid: # Reduced Gaussian grid. 1239 nlons = 2*nlats 1240 dlon = 360./nlons 1241 else: 1242 nlons = self.nx 1243 dlon = self.gridlengthXDirection 1244 lons = np.linspace(lon1,lon2,nlons) 1245 # Compute Gaussian lats (north to south) 1246 lats = gaussian_latitudes(nlats) 1247 if lat1 < lat2: # reverse them if necessary 1248 lats = lats[::-1] 1249 lons,lats = np.meshgrid(lons,lats) 1250 elif gdtn in {10,20,30,31,110}: 1251 # Mercator, Lambert Conformal, Stereographic, Albers Equal Area, 1252 # Azimuthal Equidistant 1253 dx,dy = self.gridlengthXDirection, self.gridlengthYDirection 1254 lon1,lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1255 pj = pyproj.Proj(self.projParameters) 1256 llcrnrx, llcrnry = pj(lon1,lat1) 1257 x = llcrnrx+dx*np.arange(self.nx) 1258 y = llcrnry+dy*np.arange(self.ny) 1259 x,y = np.meshgrid(x, y) 1260 lons,lats = pj(x, y, inverse=True) 1261 elif gdtn == 90: 1262 # Satellite Projection 1263 dx = self.gridlengthXDirection 1264 dy = self.gridlengthYDirection 1265 pj = pyproj.Proj(self.projParameters) 1266 x = dx*np.indices((self.ny,self.nx),'f')[1,:,:] 1267 x -= 0.5*x.max() 1268 y = dy*np.indices((self.ny,self.nx),'f')[0,:,:] 1269 y -= 0.5*y.max() 1270 lons,lats = pj(x,y,inverse=True) 1271 # Set lons,lats to 1.e30 where undefined 1272 abslons = np.fabs(lons) 1273 abslats = np.fabs(lats) 1274 lons = np.where(abslons < 1.e20, lons, 1.e30) 1275 lats = np.where(abslats < 1.e20, lats, 1.e30) 1276 elif gdtn == 32769: 1277 # Special NCEP Grid, Rotated Lat/Lon, Arakawa E-Grid (Non-Staggered) 1278 from grib2io.utils import arakawa_rotated_grid 1279 from grib2io.utils.rotated_grid import DEG2RAD 1280 di, dj = 0.0, 0.0 1281 do_180 = False 1282 idir = 1 if self.scanModeFlags[0] == 0 else -1 1283 jdir = -1 if self.scanModeFlags[1] == 0 else 1 1284 grid_oriented = 0 if self.resolutionAndComponentFlags[4] == 0 else 1 1285 do_rot = 1 if grid_oriented == 1 else 0 1286 la1 = self.latitudeFirstGridpoint 1287 lo1 = self.longitudeFirstGridpoint 1288 clon = self.longitudeCenterGridpoint 1289 clat = self.latitudeCenterGridpoint 1290 lasp = clat - 90.0 1291 losp = clon 1292 llat, llon = arakawa_rotated_grid.ll2rot(la1,lo1,lasp,losp) 1293 la2, lo2 = arakawa_rotated_grid.rot2ll(-llat,-llon,lasp,losp) 1294 rlat = -llat 1295 rlon = -llon 1296 if self.nx == 1: 1297 di = 0.0 1298 elif idir == 1: 1299 ti = rlon 1300 while ti < llon: 1301 ti += 360.0 1302 di = (ti - llon)/float(self.nx-1) 1303 else: 1304 ti = llon 1305 while ti < rlon: 1306 ti += 360.0 1307 di = (ti - rlon)/float(self.nx-1) 1308 if self.ny == 1: 1309 dj = 0.0 1310 else: 1311 dj = (rlat - llat)/float(self.ny-1) 1312 if dj < 0.0: 1313 dj = -dj 1314 if idir == 1: 1315 if llon > rlon: 1316 llon -= 360.0 1317 if llon < 0 and rlon > 0: 1318 do_180 = True 1319 else: 1320 if rlon > llon: 1321 rlon -= 360.0 1322 if rlon < 0 and llon > 0: 1323 do_180 = True 1324 xlat1d = llat + (np.arange(self.ny)*jdir*dj) 1325 xlon1d = llon + (np.arange(self.nx)*idir*di) 1326 xlons, xlats = np.meshgrid(xlon1d,xlat1d) 1327 rot2ll_vectorized = np.vectorize(arakawa_rotated_grid.rot2ll) 1328 lats, lons = rot2ll_vectorized(xlats,xlons,lasp,losp) 1329 if do_180: 1330 lons = np.where(lons>180.0,lons-360.0,lons) 1331 vector_rotation_angles_vectorized = np.vectorize(arakawa_rotated_grid.vector_rotation_angles) 1332 rots = vector_rotation_angles_vectorized(lats, lons, clat, losp, xlats) 1333 del xlat1d, xlon1d, xlats, xlons 1334 else: 1335 raise ValueError('Unsupported grid') 1336 1337 _latlon_datastore[self._sha1_section3] = dict(latitude=lats,longitude=lons) 1338 try: 1339 _latlon_datastore[self._sha1_section3]['vector_rotation_angles'] = rots 1340 except(NameError): 1341 pass 1342 1343 return lats, lons
Return lats,lons (in degrees) of grid.
Currently can handle reg. lat/lon,cglobal Gaussian, mercator, stereographic, lambert conformal, albers equal-area, space-view and azimuthal equidistant grids.
Parameters
- unrotate: If
True
[DEFAULT], and grid is rotated lat/lon, then unrotate the grid, otherwiseFalse
, do not.
Returns
- lats, lons (numpy.ndarray): Returns two numpy.ndarrays with dtype=numpy.float32 of grid latitudes and longitudes in units of degrees.
1346 def map_keys(self): 1347 """ 1348 Unpack data grid replacing integer values with strings. 1349 1350 These types of fields are categorical or classifications where data 1351 values do not represent an observable or predictable physical quantity. 1352 An example of such a field would be [Dominant Precipitation Type - 1353 DPTYPE](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-201.shtml) 1354 1355 Returns 1356 ------- 1357 map_keys 1358 numpy.ndarray of string values per element. 1359 """ 1360 hold_auto_nans = _AUTO_NANS 1361 set_auto_nans(False) 1362 1363 if (np.all(self.section1[0:2] == [7, 14]) and self.shortName == "PWTHER") or ( 1364 self._isNDFD and self.shortName in {"WX", "WWA"} 1365 ): 1366 keys = utils.decode_wx_strings(self.section2) 1367 if hasattr(self,'priMissingValue') and self.priMissingValue not in [None,0]: 1368 keys[int(self.priMissingValue)] = 'Missing' 1369 if hasattr(self,'secMissingValue') and self.secMissingValue not in [None,0]: 1370 keys[int(self.secMissingValue)] = 'Missing' 1371 u,inv = np.unique(self.data,return_inverse=True) 1372 fld = np.array([keys[x] for x in u])[inv].reshape(self.data.shape) 1373 else: 1374 # For data whose units are defined in a code table (i.e. classification or mask) 1375 tblname = re.findall(r'\d\.\d+',self.units,re.IGNORECASE)[0] 1376 fld = self.data.astype(np.int32).astype(str) 1377 tbl = tables.get_table(tblname,expand=True) 1378 for val in np.unique(fld): 1379 fld = np.where(fld==val,tbl[val],fld) 1380 set_auto_nans(hold_auto_nans) 1381 return fld
Unpack data grid replacing integer values with strings.
These types of fields are categorical or classifications where data values do not represent an observable or predictable physical quantity. An example of such a field would be Dominant Precipitation Type - DPTYPE
Returns
- map_keys: numpy.ndarray of string values per element.
1384 def to_bytes(self, validate: bool=True): 1385 """ 1386 Return packed GRIB2 message in bytes format. 1387 1388 This will be useful for exporting data in non-file formats. For example, 1389 can be used to output grib data directly to S3 using the boto3 client 1390 without the need to write a temporary file to upload first. 1391 1392 Parameters 1393 ---------- 1394 validate: default=True 1395 If `True` (DEFAULT), validates first/last four bytes for proper 1396 formatting, else returns None. If `False`, message is output as is. 1397 1398 Returns 1399 ------- 1400 to_bytes 1401 Returns GRIB2 formatted message as bytes. 1402 """ 1403 if hasattr(self,'_msg'): 1404 if validate: 1405 if self.validate(): 1406 return self._msg 1407 else: 1408 return None 1409 else: 1410 return self._msg 1411 else: 1412 return None
Return packed GRIB2 message in bytes format.
This will be useful for exporting data in non-file formats. For example, can be used to output grib data directly to S3 using the boto3 client without the need to write a temporary file to upload first.
Parameters
- validate (default=True):
If
True
(DEFAULT), validates first/last four bytes for proper formatting, else returns None. IfFalse
, message is output as is.
Returns
- to_bytes: Returns GRIB2 formatted message as bytes.
1415 def interpolate(self, method, grid_def_out, method_options=None, drtn=None, 1416 num_threads=1): 1417 """ 1418 Grib2Message Interpolator 1419 1420 Performs spatial interpolation via the [NCEPLIBS-ip 1421 library](https://github.com/NOAA-EMC/NCEPLIBS-ip). This interpolate 1422 method only supports scalar interpolation. If you need to perform 1423 vector interpolation, use the module-level `grib2io.interpolate` 1424 function. 1425 1426 Parameters 1427 ---------- 1428 method 1429 Interpolate method to use. This can either be an integer or string 1430 using the following mapping: 1431 1432 | Interpolate Scheme | Integer Value | 1433 | :---: | :---: | 1434 | 'bilinear' | 0 | 1435 | 'bicubic' | 1 | 1436 | 'neighbor' | 2 | 1437 | 'budget' | 3 | 1438 | 'spectral' | 4 | 1439 | 'neighbor-budget' | 6 | 1440 1441 grid_def_out : grib2io.Grib2GridDef 1442 Grib2GridDef object of the output grid. 1443 method_options : list of ints, optional 1444 Interpolation options. See the NCEPLIBS-ip documentation for 1445 more information on how these are used. 1446 drtn 1447 Data Representation Template to be used for the returned 1448 interpolated GRIB2 message. When `None`, the data representation 1449 template of the source GRIB2 message is used. Once again, it is the 1450 user's responsibility to properly set the Data Representation 1451 Template attributes. 1452 num_threads : int, optional 1453 Number of OpenMP threads to use for interpolation. The default 1454 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 1455 was not built with OpenMP, then this keyword argument and value 1456 will have no impact. 1457 1458 Returns 1459 ------- 1460 interpolate 1461 If interpolating to a grid, a new Grib2Message object is returned. 1462 The GRIB2 metadata of the new Grib2Message object is identical to 1463 the input except where required to be different because of the new 1464 grid specs and possibly a new data representation template. 1465 1466 If interpolating to station points, the interpolated data values are 1467 returned as a numpy.ndarray. 1468 """ 1469 section0 = self.section0 1470 section0[-1] = 0 1471 gds = [0, grid_def_out.npoints, 0, 255, grid_def_out.gdtn] 1472 section3 = np.concatenate((gds,grid_def_out.gdt)) 1473 drtn = self.drtn if drtn is None else drtn 1474 1475 msg = Grib2Message(section0,self.section1,self.section2,section3, 1476 self.section4,None,self.bitMapFlag.value,drtn=drtn) 1477 1478 msg._msgnum = -1 1479 msg._deflist = self._deflist 1480 msg._coordlist = self._coordlist 1481 shape = (msg.ny,msg.nx) 1482 ndim = 2 1483 if msg.typeOfValues == 0: 1484 dtype = 'float32' 1485 elif msg.typeOfValues == 1: 1486 dtype = 'int32' 1487 msg._data = interpolate(self.data,method,Grib2GridDef.from_section3(self.section3),grid_def_out, 1488 method_options=method_options,num_threads=num_threads).reshape(msg.ny,msg.nx) 1489 msg.section5[0] = grid_def_out.npoints 1490 return msg
Grib2Message Interpolator
Performs spatial interpolation via the NCEPLIBS-ip
library. This interpolate
method only supports scalar interpolation. If you need to perform
vector interpolation, use the module-level grib2io.interpolate
function.
Parameters
- method: Interpolate method to use. This can either be an integer or string using the following mapping:
Interpolate Scheme | Integer Value |
---|---|
'bilinear' | 0 |
'bicubic' | 1 |
'neighbor' | 2 |
'budget' | 3 |
'spectral' | 4 |
'neighbor-budget' | 6 |
- grid_def_out (grib2io.Grib2GridDef): Grib2GridDef object of the output grid.
- method_options (list of ints, optional): Interpolation options. See the NCEPLIBS-ip documentation for more information on how these are used.
- drtn: Data Representation Template to be used for the returned
interpolated GRIB2 message. When
None
, the data representation template of the source GRIB2 message is used. Once again, it is the user's responsibility to properly set the Data Representation Template attributes. - num_threads (int, optional): Number of OpenMP threads to use for interpolation. The default value is 1. If NCEPLIBS-ip and grib2io's iplib extension module was not built with OpenMP, then this keyword argument and value will have no impact.
Returns
- interpolate: If interpolating to a grid, a new Grib2Message object is returned. The GRIB2 metadata of the new Grib2Message object is identical to the input except where required to be different because of the new grid specs and possibly a new data representation template.
If interpolating to station points, the interpolated data values are returned as a numpy.ndarray.
1493 def subset(self, lats, lons): 1494 """ 1495 Return a spatial subset. 1496 1497 Currently only supports regular grids of the following types: 1498 1499 | Grid Type | gdtn | 1500 | :---: | :---: | 1501 | Latitude/Longitude, Equidistant Cylindrical, or Plate Carree | 0 | 1502 | Rotated Latitude/Longitude | 1 | 1503 | Mercator | 10 | 1504 | Polar Stereographic | 20 | 1505 | Lambert Conformal | 30 | 1506 | Albers Equal-Area | 31 | 1507 | Gaussian Latitude/Longitude | 40 | 1508 | Equatorial Azimuthal Equidistant Projection | 110 | 1509 1510 Parameters 1511 ---------- 1512 lats 1513 List or tuple of latitudes. The minimum and maximum latitudes will 1514 be used to define the southern and northern boundaries. 1515 1516 The order of the latitudes is not important. The function will 1517 determine which is the minimum and maximum. 1518 1519 The latitudes should be in decimal degrees with 0.0 at the equator, 1520 positive values in the northern hemisphere increasing to 90, and 1521 negative values in the southern hemisphere decreasing to -90. 1522 lons 1523 List or tuple of longitudes. The minimum and maximum longitudes 1524 will be used to define the western and eastern boundaries. 1525 1526 The order of the longitudes is not important. The function will 1527 determine which is the minimum and maximum. 1528 1529 GRIB2 longitudes should be in decimal degrees with 0.0 at the prime 1530 meridian, positive values increasing eastward to 360. There are no 1531 negative GRIB2 longitudes. 1532 1533 The typical west longitudes that start at 0.0 at the prime meridian 1534 and decrease to -180 westward, are converted to GRIB2 longitudes by 1535 '360 - (absolute value of the west longitude)' where typical 1536 eastern longitudes are unchanged as GRIB2 longitudes. 1537 1538 Returns 1539 ------- 1540 subset 1541 A spatial subset of a GRIB2 message. 1542 """ 1543 if self.gdtn not in [0, 1, 10, 20, 30, 31, 40, 110]: 1544 raise ValueError( 1545 """ 1546 1547Subset only works for 1548 Latitude/Longitude, Equidistant Cylindrical, or Plate Carree (gdtn=0) 1549 Rotated Latitude/Longitude (gdtn=1) 1550 Mercator (gdtn=10) 1551 Polar Stereographic (gdtn=20) 1552 Lambert Conformal (gdtn=30) 1553 Albers Equal-Area (gdtn=31) 1554 Gaussian Latitude/Longitude (gdtn=40) 1555 Equatorial Azimuthal Equidistant Projection (gdtn=110) 1556 1557""" 1558 ) 1559 1560 if self.nx == 0 or self.ny == 0: 1561 raise ValueError( 1562 """ 1563 1564Subset only works for regular grids. 1565 1566""" 1567 ) 1568 1569 newmsg = Grib2Message( 1570 np.copy(self.section0), 1571 np.copy(self.section1), 1572 np.copy(self.section2), 1573 np.copy(self.section3), 1574 np.copy(self.section4), 1575 np.copy(self.section5), 1576 ) 1577 1578 msglats, msglons = self.grid() 1579 1580 la1 = np.max(lats) 1581 lo1 = np.min(lons) 1582 la2 = np.min(lats) 1583 lo2 = np.max(lons) 1584 1585 # Find the indices of the first and last grid points to the nearest 1586 # lat/lon values in the grid. 1587 first_lat = np.abs(msglats - la1) 1588 first_lon = np.abs(msglons - lo1) 1589 max_idx = np.maximum(first_lat, first_lon) 1590 first_j, first_i = np.where(max_idx == np.min(max_idx)) 1591 1592 last_lat = np.abs(msglats - la2) 1593 last_lon = np.abs(msglons - lo2) 1594 max_idx = np.maximum(last_lat, last_lon) 1595 last_j, last_i = np.where(max_idx == np.min(max_idx)) 1596 1597 setattr(newmsg, "latitudeFirstGridpoint", msglats[first_j[0], first_i[0]]) 1598 setattr(newmsg, "longitudeFirstGridpoint", msglons[first_j[0], first_i[0]]) 1599 setattr(newmsg, "nx", np.abs(first_i[0] - last_i[0])) 1600 setattr(newmsg, "ny", np.abs(first_j[0] - last_j[0])) 1601 1602 # Set *LastGridpoint attributes even if only used for gdtn=[0, 1, 40]. 1603 # This information is used to subset xarray datasets and even though 1604 # unnecessary for some supported grid types, it won't affect a grib2io 1605 # message to set them. 1606 setattr(newmsg, "latitudeLastGridpoint", msglats[last_j[0], last_i[0]]) 1607 setattr(newmsg, "longitudeLastGridpoint", msglons[last_j[0], last_i[0]]) 1608 1609 setattr( 1610 newmsg, 1611 "data", 1612 self.data[ 1613 min(first_j[0], last_j[0]) : max(first_j[0], last_j[0]), 1614 min(first_i[0], last_i[0]) : max(first_i[0], last_i[0]), 1615 ].copy(), 1616 ) 1617 1618 # Need to set the newmsg._sha1_section3 to a blank string so the grid 1619 # method ignores the cached lat/lon values. This will force the grid 1620 # method to recompute the lat/lon values for the subsetted grid. 1621 newmsg._sha1_section3 = "" 1622 newmsg.grid() 1623 1624 return newmsg
Return a spatial subset.
Currently only supports regular grids of the following types:
Grid Type | gdtn |
---|---|
Latitude/Longitude, Equidistant Cylindrical, or Plate Carree | 0 |
Rotated Latitude/Longitude | 1 |
Mercator | 10 |
Polar Stereographic | 20 |
Lambert Conformal | 30 |
Albers Equal-Area | 31 |
Gaussian Latitude/Longitude | 40 |
Equatorial Azimuthal Equidistant Projection | 110 |
Parameters
- lats: List or tuple of latitudes. The minimum and maximum latitudes will be used to define the southern and northern boundaries.
The order of the latitudes is not important. The function will determine which is the minimum and maximum.
The latitudes should be in decimal degrees with 0.0 at the equator, positive values in the northern hemisphere increasing to 90, and negative values in the southern hemisphere decreasing to -90.
- lons: List or tuple of longitudes. The minimum and maximum longitudes will be used to define the western and eastern boundaries.
The order of the longitudes is not important. The function will determine which is the minimum and maximum.
GRIB2 longitudes should be in decimal degrees with 0.0 at the prime meridian, positive values increasing eastward to 360. There are no negative GRIB2 longitudes.
The typical west longitudes that start at 0.0 at the prime meridian and decrease to -180 westward, are converted to GRIB2 longitudes by '360 - (absolute value of the west longitude)' where typical eastern longitudes are unchanged as GRIB2 longitudes.
Returns
- subset: A spatial subset of a GRIB2 message.
1626 def validate(self): 1627 """ 1628 Validate a complete GRIB2 message. 1629 1630 The g2c library does its own internal validation when g2_gribend() is called, but 1631 we will check in grib2io also. The validation checks if the first 4 bytes in 1632 self._msg is 'GRIB' and '7777' as the last 4 bytes and that the message length in 1633 section 0 equals the length of the packed message. 1634 1635 Returns 1636 ------- 1637 `True` if the packed GRIB2 message is complete and well-formed, `False` otherwise. 1638 """ 1639 valid = False 1640 if hasattr(self,'_msg'): 1641 if self._msg[0:4]+self._msg[-4:] == b'GRIB7777': 1642 if self.section0[-1] == len(self._msg): 1643 valid = True 1644 return valid
Validate a complete GRIB2 message.
The g2c library does its own internal validation when g2_gribend() is called, but we will check in grib2io also. The validation checks if the first 4 bytes in self._msg is 'GRIB' and '7777' as the last 4 bytes and that the message length in section 0 equals the length of the packed message.
Returns
True
if the packed GRIB2 message is complete and well-formed,False
otherwise.
2119@dataclass 2120class Grib2GridDef: 2121 """ 2122 Class for Grid Definition Template Number and Template as attributes. 2123 2124 This allows for cleaner looking code when passing these metadata around. 2125 For example, the `grib2io._Grib2Message.interpolate` method and 2126 `grib2io.interpolate` function accepts these objects. 2127 """ 2128 gdtn: int 2129 gdt: NDArray 2130 2131 @classmethod 2132 def from_section3(cls, section3): 2133 return cls(section3[4],section3[5:]) 2134 2135 @property 2136 def nx(self): 2137 """Number of grid points in x-direction.""" 2138 return int(self.gdt[7]) 2139 2140 @property 2141 def ny(self): 2142 """Number of grid points in y-direction.""" 2143 return int(self.gdt[8]) 2144 2145 @property 2146 def npoints(self): 2147 """Total number of grid points.""" 2148 return int(self.gdt[7] * self.gdt[8]) 2149 2150 @property 2151 def shape(self): 2152 """Shape of the grid.""" 2153 return (int(self.ny), int(self.nx)) 2154 2155 def to_section3(self): 2156 """Return a full GRIB2 section3 array.""" 2157 return np.array( 2158 [0, self.npoints, 0, 0, self.gdtn] + list(self.gdt) 2159 ).astype(np.int64)
Class for Grid Definition Template Number and Template as attributes.
This allows for cleaner looking code when passing these metadata around.
For example, the grib2io._Grib2Message.interpolate
method and
grib2io.interpolate
function accepts these objects.
2135 @property 2136 def nx(self): 2137 """Number of grid points in x-direction.""" 2138 return int(self.gdt[7])
Number of grid points in x-direction.
2140 @property 2141 def ny(self): 2142 """Number of grid points in y-direction.""" 2143 return int(self.gdt[8])
Number of grid points in y-direction.
2145 @property 2146 def npoints(self): 2147 """Total number of grid points.""" 2148 return int(self.gdt[7] * self.gdt[8])
Total number of grid points.
653def msgs_from_index(index, filehandle=None): 654 """ 655 return list of Grib2Messages from index 656 msgs can get the data so long as an open filehandle is provided 657 """ 658 659 zipped = zip( 660 index["section0"], 661 index["section1"], 662 index["section2"], 663 index["section3"], 664 index["section4"], 665 index["section5"], 666 index["bmapflag"] 667 ) 668 msgs = [Grib2Message(*sections) for sections in zipped] 669 670 if filehandle is not None: 671 for n, (msg, offset, secpos) in enumerate(zip(msgs, index["offset"], index["sectionOffset"])): 672 msg._ondiskarray = Grib2MessageOnDiskArray(shape=(msg.ny,msg.nx), 673 ndim=2, 674 dtype=TYPE_OF_VALUES_DTYPE[msg.typeOfValues], 675 filehandle=filehandle, 676 msg=msg, 677 offset=offset, 678 bitmap_offset=secpos[6], 679 data_offset=secpos[7]) 680 msg._msgnum = n 681 return msgs
return list of Grib2Messages from index msgs can get the data so long as an open filehandle is provided