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.
1828def interpolate(a, 1829 method: Union[int, str], 1830 grid_def_in, 1831 grid_def_out, 1832 method_options=None, 1833 num_threads=1, 1834): 1835 """ 1836 This is the module-level interpolation function. 1837 1838 This interfaces with the 4-byte (32-bit float) [NCEPLIBS-ip library](https://github.com/NOAA-EMC/NCEPLIBS-ip) 1839 through grib2io's internal iplib Cython extension module. 1840 1841 Parameters 1842 ---------- 1843 a : numpy.ndarray or tuple 1844 Input data. If `a` is a `numpy.ndarray`, scalar interpolation will be 1845 performed. If `a` is a `tuple`, then vector interpolation will be 1846 performed with the assumption that u = a[0] and v = a[1] and are both 1847 `numpy.ndarray`. 1848 1849 These data are expected to be in 2-dimensional form with shape (ny, nx) 1850 or 3-dimensional (:, ny, nx) where the 1st dimension represents another 1851 spatial, temporal, or classification (i.e. ensemble members) dimension. 1852 The function will properly flatten the (ny,nx) dimensions into (nx * ny) 1853 acceptable for input into the interpolation subroutines. If needed, these 1854 data will be converted to `np.float32`. 1855 method 1856 Interpolate method to use. This can either be an integer or string using 1857 the following mapping: 1858 1859 | Interpolate Scheme | Integer Value | 1860 | :---: | :---: | 1861 | 'bilinear' | 0 | 1862 | 'bicubic' | 1 | 1863 | 'neighbor' | 2 | 1864 | 'budget' | 3 | 1865 | 'spectral' | 4 | 1866 | 'neighbor-budget' | 6 | 1867 1868 grid_def_in : grib2io.Grib2GridDef 1869 Grib2GridDef object for the input grid. 1870 grid_def_out : grib2io.Grib2GridDef 1871 Grib2GridDef object for the output grid or station points. 1872 method_options : list of ints, optional 1873 Interpolation options. See the NCEPLIBS-ip documentation for 1874 more information on how these are used. 1875 num_threads : int, optional 1876 Number of OpenMP threads to use for interpolation. The default 1877 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 1878 was not built with OpenMP, then this keyword argument and value 1879 will have no impact. 1880 1881 Returns 1882 ------- 1883 interpolate 1884 Returns a `numpy.ndarray` of dtype `np.float32` when scalar interpolation 1885 is performed or a `tuple` of `numpy.ndarray`s when vector interpolation is 1886 performed with the assumptions that 0-index is the interpolated u and 1887 1-index is the interpolated v. 1888 """ 1889 1890 from . import iplib 1891 1892 prev_num_threads = 1 1893 try: 1894 prev_num_threads = iplib.openmp_get_num_threads() 1895 iplib.openmp_set_num_threads(num_threads) 1896 except(AttributeError): 1897 pass 1898 1899 print(f"grib2io.interpolate thread report: OpenMP num threads = {iplib.openmp_get_num_threads()}") 1900 1901 if isinstance(method,int) and method not in _interp_schemes.values(): 1902 raise ValueError('Invalid interpolation method.') 1903 elif isinstance(method,str): 1904 if method in _interp_schemes.keys(): 1905 method = _interp_schemes[method] 1906 else: 1907 raise ValueError('Invalid interpolation method.') 1908 1909 if method_options is None: 1910 method_options = np.zeros((20),dtype=np.int32) 1911 if method in {3,6}: 1912 method_options[0:2] = -1 1913 1914 mi = grid_def_in.npoints 1915 mo = grid_def_out.npoints 1916 1917 # Adjust shape of input array(s) 1918 a, newshp = _adjust_array_shape_for_interp(a, grid_def_in, grid_def_out) 1919 1920 # Call interpolation subroutines according to type of a. 1921 if isinstance(a,np.ndarray): 1922 # Scalar 1923 km = a.shape[0] 1924 if np.any(np.isnan(a)): 1925 ibi = np.ones((km), dtype=np.int32) 1926 li = np.where(np.isnan(a),0,1).astype(np.uint8) 1927 else: 1928 ibi = np.zeros((km), dtype=np.int32) 1929 li = np.zeros(a.shape,dtype=np.uint8) 1930 no, rlat, rlon, ibo, lo, go, iret = iplib.interpolate_scalar( 1931 method, 1932 method_options, 1933 grid_def_in.gdtn, 1934 np.array(grid_def_in.gdt, dtype=np.int32), 1935 grid_def_out.gdtn, 1936 np.array(grid_def_out.gdt, dtype=np.int32), 1937 mi, 1938 mo, 1939 km, 1940 ibi, 1941 li, 1942 a, 1943 ) 1944 out = np.where(lo==0,np.nan,go).reshape(newshp) 1945 elif isinstance(a,tuple): 1946 # Vector 1947 km = a[0].shape[0] 1948 if np.any(np.isnan(a)): 1949 ibi = np.ones((km), dtype=np.int32) 1950 li = np.where(np.isnan(a),0,1).astype(np.uint8) 1951 else: 1952 ibi = np.zeros((km), dtype=np.int32) 1953 li = np.zeros(a[0].shape,dtype=np.uint8) 1954 no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret = iplib.interpolate_vector( 1955 method, 1956 method_options, 1957 grid_def_in.gdtn, 1958 np.array(grid_def_in.gdt, dtype=np.int32), 1959 grid_def_out.gdtn, 1960 np.array(grid_def_out.gdt, dtype=np.int32), 1961 mi, 1962 mo, 1963 km, 1964 ibi, 1965 li, 1966 a[0], 1967 a[1], 1968 ) 1969 uo = np.where(lo==0,np.nan,uo).reshape(newshp) 1970 vo = np.where(lo==0,np.nan,vo).reshape(newshp) 1971 out = (uo,vo) 1972 1973 try: 1974 iplib.openmp_set_num_threads(prev_num_threads) 1975 except(AttributeError): 1976 pass 1977 1978 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.
1981def interpolate_to_stations( 1982 a, 1983 method: Union[int, str], 1984 grid_def_in, 1985 lats, 1986 lons, 1987 method_options=None, 1988 num_threads=1 1989): 1990 """ 1991 Module-level interpolation function for interpolation to stations. 1992 1993 Interfaces with the 4-byte (32-bit float) [NCEPLIBS-ip library](https://github.com/NOAA-EMC/NCEPLIBS-ip) 1994 via grib2io's iplib Cython exntension module. It supports scalar and 1995 vector interpolation according to the type of object a. 1996 1997 Parameters 1998 ---------- 1999 a : numpy.ndarray or tuple 2000 Input data. If `a` is a `numpy.ndarray`, scalar interpolation will be 2001 performed. If `a` is a `tuple`, then vector interpolation will be 2002 performed with the assumption that u = a[0] and v = a[1] and are both 2003 `numpy.ndarray`. 2004 2005 These data are expected to be in 2-dimensional form with shape (ny, nx) 2006 or 3-dimensional (:, ny, nx) where the 1st dimension represents another 2007 spatial, temporal, or classification (i.e. ensemble members) dimension. 2008 The function will properly flatten the (ny,nx) dimensions into (nx * ny) 2009 acceptable for input into the interpolation subroutines. If needed, these 2010 data will be converted to `np.float32`. 2011 method 2012 Interpolate method to use. This can either be an integer or string using 2013 the following mapping: 2014 2015 | Interpolate Scheme | Integer Value | 2016 | :---: | :---: | 2017 | 'bilinear' | 0 | 2018 | 'bicubic' | 1 | 2019 | 'neighbor' | 2 | 2020 | 'budget' | 3 | 2021 | 'spectral' | 4 | 2022 | 'neighbor-budget' | 6 | 2023 2024 grid_def_in : grib2io.Grib2GridDef 2025 Grib2GridDef object for the input grid. 2026 lats : numpy.ndarray or list 2027 Latitudes for station points 2028 lons : numpy.ndarray or list 2029 Longitudes for station points 2030 method_options : list of ints, optional 2031 Interpolation options. See the NCEPLIBS-ip documentation for 2032 more information on how these are used. 2033 num_threads : int, optional 2034 Number of OpenMP threads to use for interpolation. The default 2035 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 2036 was not built with OpenMP, then this keyword argument and value 2037 will have no impact. 2038 2039 Returns 2040 ------- 2041 interpolate_to_stations 2042 Returns a `numpy.ndarray` of dtype `np.float32` when scalar 2043 interpolation is performed or a `tuple` of `numpy.ndarray`s 2044 when vector interpolation is performed with the assumptions 2045 that 0-index is the interpolated u and 1-index is the 2046 interpolated v. 2047 """ 2048 from . import iplib 2049 2050 prev_num_threads = 1 2051 try: 2052 prev_num_threads = iplib.openmp_get_num_threads() 2053 iplib.openmp_set_num_threads(num_threads) 2054 except(AttributeError): 2055 pass 2056 2057 if isinstance(method,int) and method not in _interp_schemes.values(): 2058 raise ValueError('Invalid interpolation method.') 2059 elif isinstance(method,str): 2060 if method in _interp_schemes.keys(): 2061 method = _interp_schemes[method] 2062 else: 2063 raise ValueError('Invalid interpolation method.') 2064 2065 if method_options is None: 2066 method_options = np.zeros((20),dtype=np.int32) 2067 if method in {3,6}: 2068 method_options[0:2] = -1 2069 2070 # Check lats and lons 2071 if isinstance(lats,list): 2072 nlats = len(lats) 2073 elif isinstance(lats,np.ndarray) and len(lats.shape) == 1: 2074 nlats = lats.shape[0] 2075 else: 2076 raise ValueError('Station latitudes must be a list or 1-D NumPy array.') 2077 if isinstance(lons,list): 2078 nlons = len(lons) 2079 elif isinstance(lons,np.ndarray) and len(lons.shape) == 1: 2080 nlons = lons.shape[0] 2081 else: 2082 raise ValueError('Station longitudes must be a list or 1-D NumPy array.') 2083 if nlats != nlons: 2084 raise ValueError('Station lats and lons must be same size.') 2085 2086 mi = grid_def_in.npoints 2087 mo = nlats 2088 2089 # Adjust shape of input array(s) 2090 a, newshp = _adjust_array_shape_for_interp_stations(a, grid_def_in, mo) 2091 2092 # Use gdtn = -1 for stations and an empty template array 2093 gdtn = -1 2094 gdt = np.zeros((200),dtype=np.int32) 2095 2096 # Call interpolation subroutines according to type of a. 2097 if isinstance(a,np.ndarray): 2098 # Scalar 2099 km = a.shape[0] 2100 ibi = np.zeros((km), dtype=np.int32) 2101 li = np.zeros(a.shape,dtype=np.uint8) 2102 no, rlat, rlon, ibo, lo, go, iret = iplib.interpolate_scalar( 2103 method, 2104 method_options, 2105 grid_def_in.gdtn, 2106 np.array(grid_def_in.gdt, dtype=np.int32), 2107 gdtn, 2108 gdt, 2109 mi, 2110 mo, 2111 km, 2112 ibi, 2113 li, 2114 a, 2115 lats=np.array(lats, dtype=np.float32), 2116 lons=np.array(lons, dtype=np.float32), 2117 ) 2118 out = go.reshape(newshp) 2119 2120 elif isinstance(a,tuple): 2121 # Vector 2122 km = a[0].shape[0] 2123 ibi = np.zeros((km), dtype=np.int32) 2124 li = np.zeros(a[0].shape,dtype=np.uint8) 2125 no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret = iplib.interpolate_vector( 2126 method, 2127 method_options, 2128 grid_def_in.gdtn, 2129 np.array(grid_def_in.gdt, dtype=np.int32), 2130 gdtn, 2131 gdt, 2132 mi, 2133 mo, 2134 km, 2135 ibi, 2136 li, 2137 a[0], 2138 a[1], 2139 lats=np.array(lats, dtype=np.float32), 2140 lons=np.array(lons, dtype=np.float32), 2141 ) 2142 out = (uo.reshape(newshp),vo.reshape(newshp)) 2143 2144 try: 2145 iplib.openmp_set_num_threads(prev_num_threads) 2146 except(AttributeError): 2147 pass 2148 2149 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, arr): 1150 """ 1151 Set the internal data array, enforcing shape (ny, nx) and dtype float32. 1152 1153 If the Grid Definition Section (section 3) of Grib2Message object is 1154 not fully formed (i.e. nx, ny = 0, 0), then the shape of the data array 1155 will be used to set nx and ny of the Grib2Message object. It will be the 1156 responsibility of the user to populate the rest of the Grid Definition 1157 Section attributes. 1158 1159 Parameters 1160 ---------- 1161 arr : array_like 1162 A 2D array whose shape must match ``(self.ny, self.nx)``. 1163 It will be converted to ``float32`` and C-contiguous if needed. 1164 1165 Raises 1166 ------ 1167 ValueError 1168 If the shape of `arr` does not match the expected dimensions. 1169 """ 1170 if not isinstance(arr, np.ndarray): 1171 raise ValueError(f"Grib2Message data only supports numpy arrays") 1172 if self.nx == 0 and self.ny == 0: 1173 self.ny = arr.shape[0] 1174 self.nx = arr.shape[1] 1175 if arr.shape != (self.ny, self.nx): 1176 raise ValueError( 1177 f"Data shape mismatch: expected ({self.ny}, {self.nx}), " 1178 f"got {arr.shape}" 1179 ) 1180 # Ensure contiguous memory layout (important for C interoperability) 1181 if not arr.flags["C_CONTIGUOUS"]: 1182 arr = np.ascontiguousarray(arr, dtype=np.float32) 1183 self._data = arr 1184 1185 1186 def flush_data(self): 1187 """ 1188 Flush the unpacked data values from the Grib2Message object. 1189 1190 Notes 1191 ----- 1192 If the Grib2Message object was constructed from "scratch" (i.e. 1193 not read from file), this method will remove the data array from 1194 the object and it cannot be recovered. 1195 """ 1196 self._data = None 1197 self.bitmap = None 1198 1199 1200 def __getitem__(self, item): 1201 return self.data[item] 1202 1203 1204 def __setitem__(self, item): 1205 raise NotImplementedError('assignment of data not supported via setitem') 1206 1207 1208 def latlons(self, *args, **kwrgs): 1209 """Alias for `grib2io.Grib2Message.grid` method.""" 1210 return self.grid(*args, **kwrgs) 1211 1212 1213 def grid(self, unrotate: bool=True): 1214 """ 1215 Return lats,lons (in degrees) of grid. 1216 1217 Currently can handle reg. lat/lon,cglobal Gaussian, mercator, 1218 stereographic, lambert conformal, albers equal-area, space-view and 1219 azimuthal equidistant grids. 1220 1221 Parameters 1222 ---------- 1223 unrotate 1224 If `True` [DEFAULT], and grid is rotated lat/lon, then unrotate the 1225 grid, otherwise `False`, do not. 1226 1227 Returns 1228 ------- 1229 lats, lons : numpy.ndarray 1230 Returns two numpy.ndarrays with dtype=numpy.float32 of grid 1231 latitudes and longitudes in units of degrees. 1232 """ 1233 if self._sha1_section3 in _latlon_datastore.keys(): 1234 return (_latlon_datastore[self._sha1_section3]['latitude'], 1235 _latlon_datastore[self._sha1_section3]['longitude']) 1236 gdtn = self.gridDefinitionTemplateNumber.value 1237 gdtmpl = self.gridDefinitionTemplate 1238 reggrid = self.gridDefinitionSection[2] == 0 # This means regular 2-d grid 1239 if gdtn == 0: 1240 # Regular lat/lon grid 1241 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1242 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1243 dlon = self.gridlengthXDirection 1244 dlat = self.gridlengthYDirection 1245 if lon2 < lon1 and dlon < 0: lon1 = -lon1 1246 lats = np.linspace(lat1,lat2,self.ny) 1247 if reggrid: 1248 lons = np.linspace(lon1,lon2,self.nx) 1249 else: 1250 lons = np.linspace(lon1,lon2,self.ny*2) 1251 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1252 elif gdtn == 1: # Rotated Lat/Lon grid 1253 pj = pyproj.Proj(self.projParameters) 1254 lat1,lon1 = self.latitudeFirstGridpoint,self.longitudeFirstGridpoint 1255 lat2,lon2 = self.latitudeLastGridpoint,self.longitudeLastGridpoint 1256 if lon1 > 180.0: lon1 -= 360.0 1257 if lon2 > 180.0: lon2 -= 360.0 1258 lats = np.linspace(lat1,lat2,self.ny) 1259 lons = np.linspace(lon1,lon2,self.nx) 1260 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1261 if unrotate: 1262 from grib2io.utils import rotated_grid 1263 lats,lons = rotated_grid.unrotate(lats,lons,self.anglePoleRotation, 1264 self.latitudeSouthernPole, 1265 self.longitudeSouthernPole) 1266 elif gdtn == 40: # Gaussian grid (only works for global!) 1267 from grib2io.utils.gauss_grid import gaussian_latitudes 1268 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1269 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1270 nlats = self.ny 1271 if not reggrid: # Reduced Gaussian grid. 1272 nlons = 2*nlats 1273 dlon = 360./nlons 1274 else: 1275 nlons = self.nx 1276 dlon = self.gridlengthXDirection 1277 lons = np.linspace(lon1,lon2,nlons) 1278 # Compute Gaussian lats (north to south) 1279 lats = gaussian_latitudes(nlats) 1280 if lat1 < lat2: # reverse them if necessary 1281 lats = lats[::-1] 1282 lons,lats = np.meshgrid(lons,lats) 1283 elif gdtn in {10,20,30,31,110}: 1284 # Mercator, Lambert Conformal, Stereographic, Albers Equal Area, 1285 # Azimuthal Equidistant 1286 dx,dy = self.gridlengthXDirection, self.gridlengthYDirection 1287 lon1,lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1288 pj = pyproj.Proj(self.projParameters) 1289 llcrnrx, llcrnry = pj(lon1,lat1) 1290 x = llcrnrx+dx*np.arange(self.nx) 1291 y = llcrnry+dy*np.arange(self.ny) 1292 x,y = np.meshgrid(x, y) 1293 lons,lats = pj(x, y, inverse=True) 1294 elif gdtn == 90: 1295 # Satellite Projection 1296 dx = self.gridlengthXDirection 1297 dy = self.gridlengthYDirection 1298 pj = pyproj.Proj(self.projParameters) 1299 x = dx*np.indices((self.ny,self.nx),'f')[1,:,:] 1300 x -= 0.5*x.max() 1301 y = dy*np.indices((self.ny,self.nx),'f')[0,:,:] 1302 y -= 0.5*y.max() 1303 lons,lats = pj(x,y,inverse=True) 1304 # Set lons,lats to 1.e30 where undefined 1305 abslons = np.fabs(lons) 1306 abslats = np.fabs(lats) 1307 lons = np.where(abslons < 1.e20, lons, 1.e30) 1308 lats = np.where(abslats < 1.e20, lats, 1.e30) 1309 elif gdtn == 32769: 1310 # Special NCEP Grid, Rotated Lat/Lon, Arakawa E-Grid (Non-Staggered) 1311 from grib2io.utils import arakawa_rotated_grid 1312 from grib2io.utils.rotated_grid import DEG2RAD 1313 di, dj = 0.0, 0.0 1314 do_180 = False 1315 idir = 1 if self.scanModeFlags[0] == 0 else -1 1316 jdir = -1 if self.scanModeFlags[1] == 0 else 1 1317 grid_oriented = 0 if self.resolutionAndComponentFlags[4] == 0 else 1 1318 do_rot = 1 if grid_oriented == 1 else 0 1319 la1 = self.latitudeFirstGridpoint 1320 lo1 = self.longitudeFirstGridpoint 1321 clon = self.longitudeCenterGridpoint 1322 clat = self.latitudeCenterGridpoint 1323 lasp = clat - 90.0 1324 losp = clon 1325 llat, llon = arakawa_rotated_grid.ll2rot(la1,lo1,lasp,losp) 1326 la2, lo2 = arakawa_rotated_grid.rot2ll(-llat,-llon,lasp,losp) 1327 rlat = -llat 1328 rlon = -llon 1329 if self.nx == 1: 1330 di = 0.0 1331 elif idir == 1: 1332 ti = rlon 1333 while ti < llon: 1334 ti += 360.0 1335 di = (ti - llon)/float(self.nx-1) 1336 else: 1337 ti = llon 1338 while ti < rlon: 1339 ti += 360.0 1340 di = (ti - rlon)/float(self.nx-1) 1341 if self.ny == 1: 1342 dj = 0.0 1343 else: 1344 dj = (rlat - llat)/float(self.ny-1) 1345 if dj < 0.0: 1346 dj = -dj 1347 if idir == 1: 1348 if llon > rlon: 1349 llon -= 360.0 1350 if llon < 0 and rlon > 0: 1351 do_180 = True 1352 else: 1353 if rlon > llon: 1354 rlon -= 360.0 1355 if rlon < 0 and llon > 0: 1356 do_180 = True 1357 xlat1d = llat + (np.arange(self.ny)*jdir*dj) 1358 xlon1d = llon + (np.arange(self.nx)*idir*di) 1359 xlons, xlats = np.meshgrid(xlon1d,xlat1d) 1360 rot2ll_vectorized = np.vectorize(arakawa_rotated_grid.rot2ll) 1361 lats, lons = rot2ll_vectorized(xlats,xlons,lasp,losp) 1362 if do_180: 1363 lons = np.where(lons>180.0,lons-360.0,lons) 1364 vector_rotation_angles_vectorized = np.vectorize(arakawa_rotated_grid.vector_rotation_angles) 1365 rots = vector_rotation_angles_vectorized(lats, lons, clat, losp, xlats) 1366 del xlat1d, xlon1d, xlats, xlons 1367 else: 1368 raise ValueError('Unsupported grid') 1369 1370 _latlon_datastore[self._sha1_section3] = dict(latitude=lats,longitude=lons) 1371 try: 1372 _latlon_datastore[self._sha1_section3]['vector_rotation_angles'] = rots 1373 except(NameError): 1374 pass 1375 1376 return lats, lons 1377 1378 1379 def map_keys(self): 1380 """ 1381 Unpack data grid replacing integer values with strings. 1382 1383 These types of fields are categorical or classifications where data 1384 values do not represent an observable or predictable physical quantity. 1385 An example of such a field would be [Dominant Precipitation Type - 1386 DPTYPE](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-201.shtml) 1387 1388 Returns 1389 ------- 1390 map_keys 1391 numpy.ndarray of string values per element. 1392 """ 1393 hold_auto_nans = _AUTO_NANS 1394 set_auto_nans(False) 1395 1396 if (np.all(self.section1[0:2] == [7, 14]) and self.shortName == "PWTHER") or ( 1397 self._isNDFD and self.shortName in {"WX", "WWA"} 1398 ): 1399 keys = utils.decode_wx_strings(self.section2) 1400 if hasattr(self,'priMissingValue') and self.priMissingValue not in [None,0]: 1401 keys[int(self.priMissingValue)] = 'Missing' 1402 if hasattr(self,'secMissingValue') and self.secMissingValue not in [None,0]: 1403 keys[int(self.secMissingValue)] = 'Missing' 1404 u,inv = np.unique(self.data,return_inverse=True) 1405 fld = np.array([keys[x] for x in u])[inv].reshape(self.data.shape) 1406 else: 1407 # For data whose units are defined in a code table (i.e. classification or mask) 1408 tblname = re.findall(r'\d\.\d+',self.units,re.IGNORECASE)[0] 1409 fld = self.data.astype(np.int32).astype(str) 1410 tbl = tables.get_table(tblname,expand=True) 1411 for val in np.unique(fld): 1412 fld = np.where(fld==val,tbl[val],fld) 1413 set_auto_nans(hold_auto_nans) 1414 return fld 1415 1416 1417 def to_bytes(self, validate: bool=True): 1418 """ 1419 Return packed GRIB2 message in bytes format. 1420 1421 This will be useful for exporting data in non-file formats. For example, 1422 can be used to output grib data directly to S3 using the boto3 client 1423 without the need to write a temporary file to upload first. 1424 1425 Parameters 1426 ---------- 1427 validate: default=True 1428 If `True` (DEFAULT), validates first/last four bytes for proper 1429 formatting, else returns None. If `False`, message is output as is. 1430 1431 Returns 1432 ------- 1433 to_bytes 1434 Returns GRIB2 formatted message as bytes. 1435 """ 1436 if hasattr(self,'_msg'): 1437 if validate: 1438 if self.validate(): 1439 return self._msg 1440 else: 1441 return None 1442 else: 1443 return self._msg 1444 else: 1445 return None 1446 1447 1448 def interpolate(self, method, grid_def_out, method_options=None, drtn=None, 1449 num_threads=1): 1450 """ 1451 Grib2Message Interpolator 1452 1453 Performs spatial interpolation via the [NCEPLIBS-ip 1454 library](https://github.com/NOAA-EMC/NCEPLIBS-ip). This interpolate 1455 method only supports scalar interpolation. If you need to perform 1456 vector interpolation, use the module-level `grib2io.interpolate` 1457 function. 1458 1459 Parameters 1460 ---------- 1461 method 1462 Interpolate method to use. This can either be an integer or string 1463 using the following mapping: 1464 1465 | Interpolate Scheme | Integer Value | 1466 | :---: | :---: | 1467 | 'bilinear' | 0 | 1468 | 'bicubic' | 1 | 1469 | 'neighbor' | 2 | 1470 | 'budget' | 3 | 1471 | 'spectral' | 4 | 1472 | 'neighbor-budget' | 6 | 1473 1474 grid_def_out : grib2io.Grib2GridDef 1475 Grib2GridDef object of the output grid. 1476 method_options : list of ints, optional 1477 Interpolation options. See the NCEPLIBS-ip documentation for 1478 more information on how these are used. 1479 drtn 1480 Data Representation Template to be used for the returned 1481 interpolated GRIB2 message. When `None`, the data representation 1482 template of the source GRIB2 message is used. Once again, it is the 1483 user's responsibility to properly set the Data Representation 1484 Template attributes. 1485 num_threads : int, optional 1486 Number of OpenMP threads to use for interpolation. The default 1487 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 1488 was not built with OpenMP, then this keyword argument and value 1489 will have no impact. 1490 1491 Returns 1492 ------- 1493 interpolate 1494 If interpolating to a grid, a new Grib2Message object is returned. 1495 The GRIB2 metadata of the new Grib2Message object is identical to 1496 the input except where required to be different because of the new 1497 grid specs and possibly a new data representation template. 1498 1499 If interpolating to station points, the interpolated data values are 1500 returned as a numpy.ndarray. 1501 """ 1502 section0 = self.section0 1503 section0[-1] = 0 1504 gds = [0, grid_def_out.npoints, 0, 255, grid_def_out.gdtn] 1505 section3 = np.concatenate((gds,grid_def_out.gdt)) 1506 drtn = self.drtn if drtn is None else drtn 1507 1508 msg = Grib2Message(section0,self.section1,self.section2,section3, 1509 self.section4,None,self.bitMapFlag.value,drtn=drtn) 1510 1511 msg._msgnum = -1 1512 msg._deflist = self._deflist 1513 msg._coordlist = self._coordlist 1514 shape = (msg.ny,msg.nx) 1515 ndim = 2 1516 if msg.typeOfValues == 0: 1517 dtype = 'float32' 1518 elif msg.typeOfValues == 1: 1519 dtype = 'int32' 1520 msg._data = interpolate(self.data,method,Grib2GridDef.from_section3(self.section3),grid_def_out, 1521 method_options=method_options,num_threads=num_threads).reshape(msg.ny,msg.nx) 1522 msg.section5[0] = grid_def_out.npoints 1523 return msg 1524 1525 1526 def subset(self, lats, lons): 1527 """ 1528 Return a spatial subset. 1529 1530 Currently only supports regular grids of the following types: 1531 1532 | Grid Type | gdtn | 1533 | :---: | :---: | 1534 | Latitude/Longitude, Equidistant Cylindrical, or Plate Carree | 0 | 1535 | Rotated Latitude/Longitude | 1 | 1536 | Mercator | 10 | 1537 | Polar Stereographic | 20 | 1538 | Lambert Conformal | 30 | 1539 | Albers Equal-Area | 31 | 1540 | Gaussian Latitude/Longitude | 40 | 1541 | Equatorial Azimuthal Equidistant Projection | 110 | 1542 1543 Parameters 1544 ---------- 1545 lats 1546 List or tuple of latitudes. The minimum and maximum latitudes will 1547 be used to define the southern and northern boundaries. 1548 1549 The order of the latitudes is not important. The function will 1550 determine which is the minimum and maximum. 1551 1552 The latitudes should be in decimal degrees with 0.0 at the equator, 1553 positive values in the northern hemisphere increasing to 90, and 1554 negative values in the southern hemisphere decreasing to -90. 1555 lons 1556 List or tuple of longitudes. The minimum and maximum longitudes 1557 will be used to define the western and eastern boundaries. 1558 1559 The order of the longitudes is not important. The function will 1560 determine which is the minimum and maximum. 1561 1562 GRIB2 longitudes should be in decimal degrees with 0.0 at the prime 1563 meridian, positive values increasing eastward to 360. There are no 1564 negative GRIB2 longitudes. 1565 1566 The typical west longitudes that start at 0.0 at the prime meridian 1567 and decrease to -180 westward, are converted to GRIB2 longitudes by 1568 '360 - (absolute value of the west longitude)' where typical 1569 eastern longitudes are unchanged as GRIB2 longitudes. 1570 1571 Returns 1572 ------- 1573 subset 1574 A spatial subset of a GRIB2 message. 1575 """ 1576 if self.gdtn not in [0, 1, 10, 20, 30, 31, 40, 110]: 1577 raise ValueError( 1578 """ 1579 1580Subset only works for 1581 Latitude/Longitude, Equidistant Cylindrical, or Plate Carree (gdtn=0) 1582 Rotated Latitude/Longitude (gdtn=1) 1583 Mercator (gdtn=10) 1584 Polar Stereographic (gdtn=20) 1585 Lambert Conformal (gdtn=30) 1586 Albers Equal-Area (gdtn=31) 1587 Gaussian Latitude/Longitude (gdtn=40) 1588 Equatorial Azimuthal Equidistant Projection (gdtn=110) 1589 1590""" 1591 ) 1592 1593 if self.nx == 0 or self.ny == 0: 1594 raise ValueError( 1595 """ 1596 1597Subset only works for regular grids. 1598 1599""" 1600 ) 1601 1602 newmsg = Grib2Message( 1603 np.copy(self.section0), 1604 np.copy(self.section1), 1605 np.copy(self.section2), 1606 np.copy(self.section3), 1607 np.copy(self.section4), 1608 np.copy(self.section5), 1609 ) 1610 1611 msglats, msglons = self.grid() 1612 1613 la1 = np.max(lats) 1614 lo1 = np.min(lons) 1615 la2 = np.min(lats) 1616 lo2 = np.max(lons) 1617 1618 # Find the indices of the first and last grid points to the nearest 1619 # lat/lon values in the grid. 1620 first_lat = np.abs(msglats - la1) 1621 first_lon = np.abs(msglons - lo1) 1622 max_idx = np.maximum(first_lat, first_lon) 1623 first_j, first_i = np.where(max_idx == np.min(max_idx)) 1624 1625 last_lat = np.abs(msglats - la2) 1626 last_lon = np.abs(msglons - lo2) 1627 max_idx = np.maximum(last_lat, last_lon) 1628 last_j, last_i = np.where(max_idx == np.min(max_idx)) 1629 1630 setattr(newmsg, "latitudeFirstGridpoint", msglats[first_j[0], first_i[0]]) 1631 setattr(newmsg, "longitudeFirstGridpoint", msglons[first_j[0], first_i[0]]) 1632 setattr(newmsg, "nx", np.abs(first_i[0] - last_i[0])) 1633 setattr(newmsg, "ny", np.abs(first_j[0] - last_j[0])) 1634 1635 # Set *LastGridpoint attributes even if only used for gdtn=[0, 1, 40]. 1636 # This information is used to subset xarray datasets and even though 1637 # unnecessary for some supported grid types, it won't affect a grib2io 1638 # message to set them. 1639 setattr(newmsg, "latitudeLastGridpoint", msglats[last_j[0], last_i[0]]) 1640 setattr(newmsg, "longitudeLastGridpoint", msglons[last_j[0], last_i[0]]) 1641 1642 setattr( 1643 newmsg, 1644 "data", 1645 self.data[ 1646 min(first_j[0], last_j[0]) : max(first_j[0], last_j[0]), 1647 min(first_i[0], last_i[0]) : max(first_i[0], last_i[0]), 1648 ].copy(), 1649 ) 1650 1651 # Need to set the newmsg._sha1_section3 to a blank string so the grid 1652 # method ignores the cached lat/lon values. This will force the grid 1653 # method to recompute the lat/lon values for the subsetted grid. 1654 newmsg._sha1_section3 = "" 1655 newmsg.grid() 1656 1657 return newmsg 1658 1659 def validate(self): 1660 """ 1661 Validate a complete GRIB2 message. 1662 1663 The g2c library does its own internal validation when g2_gribend() is called, but 1664 we will check in grib2io also. The validation checks if the first 4 bytes in 1665 self._msg is 'GRIB' and '7777' as the last 4 bytes and that the message length in 1666 section 0 equals the length of the packed message. 1667 1668 Returns 1669 ------- 1670 `True` if the packed GRIB2 message is complete and well-formed, `False` otherwise. 1671 """ 1672 valid = False 1673 if hasattr(self,'_msg'): 1674 if self._msg[0:4]+self._msg[-4:] == b'GRIB7777': 1675 if self.section0[-1] == len(self._msg): 1676 valid = True 1677 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.
1186 def flush_data(self): 1187 """ 1188 Flush the unpacked data values from the Grib2Message object. 1189 1190 Notes 1191 ----- 1192 If the Grib2Message object was constructed from "scratch" (i.e. 1193 not read from file), this method will remove the data array from 1194 the object and it cannot be recovered. 1195 """ 1196 self._data = None 1197 self.bitmap = None
Flush the unpacked data values from the Grib2Message object.
Notes
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.
1208 def latlons(self, *args, **kwrgs): 1209 """Alias for `grib2io.Grib2Message.grid` method.""" 1210 return self.grid(*args, **kwrgs)
Alias for grib2io.Grib2Message.grid
method.
1213 def grid(self, unrotate: bool=True): 1214 """ 1215 Return lats,lons (in degrees) of grid. 1216 1217 Currently can handle reg. lat/lon,cglobal Gaussian, mercator, 1218 stereographic, lambert conformal, albers equal-area, space-view and 1219 azimuthal equidistant grids. 1220 1221 Parameters 1222 ---------- 1223 unrotate 1224 If `True` [DEFAULT], and grid is rotated lat/lon, then unrotate the 1225 grid, otherwise `False`, do not. 1226 1227 Returns 1228 ------- 1229 lats, lons : numpy.ndarray 1230 Returns two numpy.ndarrays with dtype=numpy.float32 of grid 1231 latitudes and longitudes in units of degrees. 1232 """ 1233 if self._sha1_section3 in _latlon_datastore.keys(): 1234 return (_latlon_datastore[self._sha1_section3]['latitude'], 1235 _latlon_datastore[self._sha1_section3]['longitude']) 1236 gdtn = self.gridDefinitionTemplateNumber.value 1237 gdtmpl = self.gridDefinitionTemplate 1238 reggrid = self.gridDefinitionSection[2] == 0 # This means regular 2-d grid 1239 if gdtn == 0: 1240 # Regular lat/lon grid 1241 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1242 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1243 dlon = self.gridlengthXDirection 1244 dlat = self.gridlengthYDirection 1245 if lon2 < lon1 and dlon < 0: lon1 = -lon1 1246 lats = np.linspace(lat1,lat2,self.ny) 1247 if reggrid: 1248 lons = np.linspace(lon1,lon2,self.nx) 1249 else: 1250 lons = np.linspace(lon1,lon2,self.ny*2) 1251 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1252 elif gdtn == 1: # Rotated Lat/Lon grid 1253 pj = pyproj.Proj(self.projParameters) 1254 lat1,lon1 = self.latitudeFirstGridpoint,self.longitudeFirstGridpoint 1255 lat2,lon2 = self.latitudeLastGridpoint,self.longitudeLastGridpoint 1256 if lon1 > 180.0: lon1 -= 360.0 1257 if lon2 > 180.0: lon2 -= 360.0 1258 lats = np.linspace(lat1,lat2,self.ny) 1259 lons = np.linspace(lon1,lon2,self.nx) 1260 lons,lats = np.meshgrid(lons,lats) # Make 2-d arrays. 1261 if unrotate: 1262 from grib2io.utils import rotated_grid 1263 lats,lons = rotated_grid.unrotate(lats,lons,self.anglePoleRotation, 1264 self.latitudeSouthernPole, 1265 self.longitudeSouthernPole) 1266 elif gdtn == 40: # Gaussian grid (only works for global!) 1267 from grib2io.utils.gauss_grid import gaussian_latitudes 1268 lon1, lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1269 lon2, lat2 = self.longitudeLastGridpoint, self.latitudeLastGridpoint 1270 nlats = self.ny 1271 if not reggrid: # Reduced Gaussian grid. 1272 nlons = 2*nlats 1273 dlon = 360./nlons 1274 else: 1275 nlons = self.nx 1276 dlon = self.gridlengthXDirection 1277 lons = np.linspace(lon1,lon2,nlons) 1278 # Compute Gaussian lats (north to south) 1279 lats = gaussian_latitudes(nlats) 1280 if lat1 < lat2: # reverse them if necessary 1281 lats = lats[::-1] 1282 lons,lats = np.meshgrid(lons,lats) 1283 elif gdtn in {10,20,30,31,110}: 1284 # Mercator, Lambert Conformal, Stereographic, Albers Equal Area, 1285 # Azimuthal Equidistant 1286 dx,dy = self.gridlengthXDirection, self.gridlengthYDirection 1287 lon1,lat1 = self.longitudeFirstGridpoint, self.latitudeFirstGridpoint 1288 pj = pyproj.Proj(self.projParameters) 1289 llcrnrx, llcrnry = pj(lon1,lat1) 1290 x = llcrnrx+dx*np.arange(self.nx) 1291 y = llcrnry+dy*np.arange(self.ny) 1292 x,y = np.meshgrid(x, y) 1293 lons,lats = pj(x, y, inverse=True) 1294 elif gdtn == 90: 1295 # Satellite Projection 1296 dx = self.gridlengthXDirection 1297 dy = self.gridlengthYDirection 1298 pj = pyproj.Proj(self.projParameters) 1299 x = dx*np.indices((self.ny,self.nx),'f')[1,:,:] 1300 x -= 0.5*x.max() 1301 y = dy*np.indices((self.ny,self.nx),'f')[0,:,:] 1302 y -= 0.5*y.max() 1303 lons,lats = pj(x,y,inverse=True) 1304 # Set lons,lats to 1.e30 where undefined 1305 abslons = np.fabs(lons) 1306 abslats = np.fabs(lats) 1307 lons = np.where(abslons < 1.e20, lons, 1.e30) 1308 lats = np.where(abslats < 1.e20, lats, 1.e30) 1309 elif gdtn == 32769: 1310 # Special NCEP Grid, Rotated Lat/Lon, Arakawa E-Grid (Non-Staggered) 1311 from grib2io.utils import arakawa_rotated_grid 1312 from grib2io.utils.rotated_grid import DEG2RAD 1313 di, dj = 0.0, 0.0 1314 do_180 = False 1315 idir = 1 if self.scanModeFlags[0] == 0 else -1 1316 jdir = -1 if self.scanModeFlags[1] == 0 else 1 1317 grid_oriented = 0 if self.resolutionAndComponentFlags[4] == 0 else 1 1318 do_rot = 1 if grid_oriented == 1 else 0 1319 la1 = self.latitudeFirstGridpoint 1320 lo1 = self.longitudeFirstGridpoint 1321 clon = self.longitudeCenterGridpoint 1322 clat = self.latitudeCenterGridpoint 1323 lasp = clat - 90.0 1324 losp = clon 1325 llat, llon = arakawa_rotated_grid.ll2rot(la1,lo1,lasp,losp) 1326 la2, lo2 = arakawa_rotated_grid.rot2ll(-llat,-llon,lasp,losp) 1327 rlat = -llat 1328 rlon = -llon 1329 if self.nx == 1: 1330 di = 0.0 1331 elif idir == 1: 1332 ti = rlon 1333 while ti < llon: 1334 ti += 360.0 1335 di = (ti - llon)/float(self.nx-1) 1336 else: 1337 ti = llon 1338 while ti < rlon: 1339 ti += 360.0 1340 di = (ti - rlon)/float(self.nx-1) 1341 if self.ny == 1: 1342 dj = 0.0 1343 else: 1344 dj = (rlat - llat)/float(self.ny-1) 1345 if dj < 0.0: 1346 dj = -dj 1347 if idir == 1: 1348 if llon > rlon: 1349 llon -= 360.0 1350 if llon < 0 and rlon > 0: 1351 do_180 = True 1352 else: 1353 if rlon > llon: 1354 rlon -= 360.0 1355 if rlon < 0 and llon > 0: 1356 do_180 = True 1357 xlat1d = llat + (np.arange(self.ny)*jdir*dj) 1358 xlon1d = llon + (np.arange(self.nx)*idir*di) 1359 xlons, xlats = np.meshgrid(xlon1d,xlat1d) 1360 rot2ll_vectorized = np.vectorize(arakawa_rotated_grid.rot2ll) 1361 lats, lons = rot2ll_vectorized(xlats,xlons,lasp,losp) 1362 if do_180: 1363 lons = np.where(lons>180.0,lons-360.0,lons) 1364 vector_rotation_angles_vectorized = np.vectorize(arakawa_rotated_grid.vector_rotation_angles) 1365 rots = vector_rotation_angles_vectorized(lats, lons, clat, losp, xlats) 1366 del xlat1d, xlon1d, xlats, xlons 1367 else: 1368 raise ValueError('Unsupported grid') 1369 1370 _latlon_datastore[self._sha1_section3] = dict(latitude=lats,longitude=lons) 1371 try: 1372 _latlon_datastore[self._sha1_section3]['vector_rotation_angles'] = rots 1373 except(NameError): 1374 pass 1375 1376 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.
1379 def map_keys(self): 1380 """ 1381 Unpack data grid replacing integer values with strings. 1382 1383 These types of fields are categorical or classifications where data 1384 values do not represent an observable or predictable physical quantity. 1385 An example of such a field would be [Dominant Precipitation Type - 1386 DPTYPE](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-201.shtml) 1387 1388 Returns 1389 ------- 1390 map_keys 1391 numpy.ndarray of string values per element. 1392 """ 1393 hold_auto_nans = _AUTO_NANS 1394 set_auto_nans(False) 1395 1396 if (np.all(self.section1[0:2] == [7, 14]) and self.shortName == "PWTHER") or ( 1397 self._isNDFD and self.shortName in {"WX", "WWA"} 1398 ): 1399 keys = utils.decode_wx_strings(self.section2) 1400 if hasattr(self,'priMissingValue') and self.priMissingValue not in [None,0]: 1401 keys[int(self.priMissingValue)] = 'Missing' 1402 if hasattr(self,'secMissingValue') and self.secMissingValue not in [None,0]: 1403 keys[int(self.secMissingValue)] = 'Missing' 1404 u,inv = np.unique(self.data,return_inverse=True) 1405 fld = np.array([keys[x] for x in u])[inv].reshape(self.data.shape) 1406 else: 1407 # For data whose units are defined in a code table (i.e. classification or mask) 1408 tblname = re.findall(r'\d\.\d+',self.units,re.IGNORECASE)[0] 1409 fld = self.data.astype(np.int32).astype(str) 1410 tbl = tables.get_table(tblname,expand=True) 1411 for val in np.unique(fld): 1412 fld = np.where(fld==val,tbl[val],fld) 1413 set_auto_nans(hold_auto_nans) 1414 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.
1417 def to_bytes(self, validate: bool=True): 1418 """ 1419 Return packed GRIB2 message in bytes format. 1420 1421 This will be useful for exporting data in non-file formats. For example, 1422 can be used to output grib data directly to S3 using the boto3 client 1423 without the need to write a temporary file to upload first. 1424 1425 Parameters 1426 ---------- 1427 validate: default=True 1428 If `True` (DEFAULT), validates first/last four bytes for proper 1429 formatting, else returns None. If `False`, message is output as is. 1430 1431 Returns 1432 ------- 1433 to_bytes 1434 Returns GRIB2 formatted message as bytes. 1435 """ 1436 if hasattr(self,'_msg'): 1437 if validate: 1438 if self.validate(): 1439 return self._msg 1440 else: 1441 return None 1442 else: 1443 return self._msg 1444 else: 1445 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.
1448 def interpolate(self, method, grid_def_out, method_options=None, drtn=None, 1449 num_threads=1): 1450 """ 1451 Grib2Message Interpolator 1452 1453 Performs spatial interpolation via the [NCEPLIBS-ip 1454 library](https://github.com/NOAA-EMC/NCEPLIBS-ip). This interpolate 1455 method only supports scalar interpolation. If you need to perform 1456 vector interpolation, use the module-level `grib2io.interpolate` 1457 function. 1458 1459 Parameters 1460 ---------- 1461 method 1462 Interpolate method to use. This can either be an integer or string 1463 using the following mapping: 1464 1465 | Interpolate Scheme | Integer Value | 1466 | :---: | :---: | 1467 | 'bilinear' | 0 | 1468 | 'bicubic' | 1 | 1469 | 'neighbor' | 2 | 1470 | 'budget' | 3 | 1471 | 'spectral' | 4 | 1472 | 'neighbor-budget' | 6 | 1473 1474 grid_def_out : grib2io.Grib2GridDef 1475 Grib2GridDef object of the output grid. 1476 method_options : list of ints, optional 1477 Interpolation options. See the NCEPLIBS-ip documentation for 1478 more information on how these are used. 1479 drtn 1480 Data Representation Template to be used for the returned 1481 interpolated GRIB2 message. When `None`, the data representation 1482 template of the source GRIB2 message is used. Once again, it is the 1483 user's responsibility to properly set the Data Representation 1484 Template attributes. 1485 num_threads : int, optional 1486 Number of OpenMP threads to use for interpolation. The default 1487 value is 1. If NCEPLIBS-ip and grib2io's iplib extension module 1488 was not built with OpenMP, then this keyword argument and value 1489 will have no impact. 1490 1491 Returns 1492 ------- 1493 interpolate 1494 If interpolating to a grid, a new Grib2Message object is returned. 1495 The GRIB2 metadata of the new Grib2Message object is identical to 1496 the input except where required to be different because of the new 1497 grid specs and possibly a new data representation template. 1498 1499 If interpolating to station points, the interpolated data values are 1500 returned as a numpy.ndarray. 1501 """ 1502 section0 = self.section0 1503 section0[-1] = 0 1504 gds = [0, grid_def_out.npoints, 0, 255, grid_def_out.gdtn] 1505 section3 = np.concatenate((gds,grid_def_out.gdt)) 1506 drtn = self.drtn if drtn is None else drtn 1507 1508 msg = Grib2Message(section0,self.section1,self.section2,section3, 1509 self.section4,None,self.bitMapFlag.value,drtn=drtn) 1510 1511 msg._msgnum = -1 1512 msg._deflist = self._deflist 1513 msg._coordlist = self._coordlist 1514 shape = (msg.ny,msg.nx) 1515 ndim = 2 1516 if msg.typeOfValues == 0: 1517 dtype = 'float32' 1518 elif msg.typeOfValues == 1: 1519 dtype = 'int32' 1520 msg._data = interpolate(self.data,method,Grib2GridDef.from_section3(self.section3),grid_def_out, 1521 method_options=method_options,num_threads=num_threads).reshape(msg.ny,msg.nx) 1522 msg.section5[0] = grid_def_out.npoints 1523 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.
1526 def subset(self, lats, lons): 1527 """ 1528 Return a spatial subset. 1529 1530 Currently only supports regular grids of the following types: 1531 1532 | Grid Type | gdtn | 1533 | :---: | :---: | 1534 | Latitude/Longitude, Equidistant Cylindrical, or Plate Carree | 0 | 1535 | Rotated Latitude/Longitude | 1 | 1536 | Mercator | 10 | 1537 | Polar Stereographic | 20 | 1538 | Lambert Conformal | 30 | 1539 | Albers Equal-Area | 31 | 1540 | Gaussian Latitude/Longitude | 40 | 1541 | Equatorial Azimuthal Equidistant Projection | 110 | 1542 1543 Parameters 1544 ---------- 1545 lats 1546 List or tuple of latitudes. The minimum and maximum latitudes will 1547 be used to define the southern and northern boundaries. 1548 1549 The order of the latitudes is not important. The function will 1550 determine which is the minimum and maximum. 1551 1552 The latitudes should be in decimal degrees with 0.0 at the equator, 1553 positive values in the northern hemisphere increasing to 90, and 1554 negative values in the southern hemisphere decreasing to -90. 1555 lons 1556 List or tuple of longitudes. The minimum and maximum longitudes 1557 will be used to define the western and eastern boundaries. 1558 1559 The order of the longitudes is not important. The function will 1560 determine which is the minimum and maximum. 1561 1562 GRIB2 longitudes should be in decimal degrees with 0.0 at the prime 1563 meridian, positive values increasing eastward to 360. There are no 1564 negative GRIB2 longitudes. 1565 1566 The typical west longitudes that start at 0.0 at the prime meridian 1567 and decrease to -180 westward, are converted to GRIB2 longitudes by 1568 '360 - (absolute value of the west longitude)' where typical 1569 eastern longitudes are unchanged as GRIB2 longitudes. 1570 1571 Returns 1572 ------- 1573 subset 1574 A spatial subset of a GRIB2 message. 1575 """ 1576 if self.gdtn not in [0, 1, 10, 20, 30, 31, 40, 110]: 1577 raise ValueError( 1578 """ 1579 1580Subset only works for 1581 Latitude/Longitude, Equidistant Cylindrical, or Plate Carree (gdtn=0) 1582 Rotated Latitude/Longitude (gdtn=1) 1583 Mercator (gdtn=10) 1584 Polar Stereographic (gdtn=20) 1585 Lambert Conformal (gdtn=30) 1586 Albers Equal-Area (gdtn=31) 1587 Gaussian Latitude/Longitude (gdtn=40) 1588 Equatorial Azimuthal Equidistant Projection (gdtn=110) 1589 1590""" 1591 ) 1592 1593 if self.nx == 0 or self.ny == 0: 1594 raise ValueError( 1595 """ 1596 1597Subset only works for regular grids. 1598 1599""" 1600 ) 1601 1602 newmsg = Grib2Message( 1603 np.copy(self.section0), 1604 np.copy(self.section1), 1605 np.copy(self.section2), 1606 np.copy(self.section3), 1607 np.copy(self.section4), 1608 np.copy(self.section5), 1609 ) 1610 1611 msglats, msglons = self.grid() 1612 1613 la1 = np.max(lats) 1614 lo1 = np.min(lons) 1615 la2 = np.min(lats) 1616 lo2 = np.max(lons) 1617 1618 # Find the indices of the first and last grid points to the nearest 1619 # lat/lon values in the grid. 1620 first_lat = np.abs(msglats - la1) 1621 first_lon = np.abs(msglons - lo1) 1622 max_idx = np.maximum(first_lat, first_lon) 1623 first_j, first_i = np.where(max_idx == np.min(max_idx)) 1624 1625 last_lat = np.abs(msglats - la2) 1626 last_lon = np.abs(msglons - lo2) 1627 max_idx = np.maximum(last_lat, last_lon) 1628 last_j, last_i = np.where(max_idx == np.min(max_idx)) 1629 1630 setattr(newmsg, "latitudeFirstGridpoint", msglats[first_j[0], first_i[0]]) 1631 setattr(newmsg, "longitudeFirstGridpoint", msglons[first_j[0], first_i[0]]) 1632 setattr(newmsg, "nx", np.abs(first_i[0] - last_i[0])) 1633 setattr(newmsg, "ny", np.abs(first_j[0] - last_j[0])) 1634 1635 # Set *LastGridpoint attributes even if only used for gdtn=[0, 1, 40]. 1636 # This information is used to subset xarray datasets and even though 1637 # unnecessary for some supported grid types, it won't affect a grib2io 1638 # message to set them. 1639 setattr(newmsg, "latitudeLastGridpoint", msglats[last_j[0], last_i[0]]) 1640 setattr(newmsg, "longitudeLastGridpoint", msglons[last_j[0], last_i[0]]) 1641 1642 setattr( 1643 newmsg, 1644 "data", 1645 self.data[ 1646 min(first_j[0], last_j[0]) : max(first_j[0], last_j[0]), 1647 min(first_i[0], last_i[0]) : max(first_i[0], last_i[0]), 1648 ].copy(), 1649 ) 1650 1651 # Need to set the newmsg._sha1_section3 to a blank string so the grid 1652 # method ignores the cached lat/lon values. This will force the grid 1653 # method to recompute the lat/lon values for the subsetted grid. 1654 newmsg._sha1_section3 = "" 1655 newmsg.grid() 1656 1657 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.
1659 def validate(self): 1660 """ 1661 Validate a complete GRIB2 message. 1662 1663 The g2c library does its own internal validation when g2_gribend() is called, but 1664 we will check in grib2io also. The validation checks if the first 4 bytes in 1665 self._msg is 'GRIB' and '7777' as the last 4 bytes and that the message length in 1666 section 0 equals the length of the packed message. 1667 1668 Returns 1669 ------- 1670 `True` if the packed GRIB2 message is complete and well-formed, `False` otherwise. 1671 """ 1672 valid = False 1673 if hasattr(self,'_msg'): 1674 if self._msg[0:4]+self._msg[-4:] == b'GRIB7777': 1675 if self.section0[-1] == len(self._msg): 1676 valid = True 1677 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.
2152@dataclass 2153class Grib2GridDef: 2154 """ 2155 Class for Grid Definition Template Number and Template as attributes. 2156 2157 This allows for cleaner looking code when passing these metadata around. 2158 For example, the `grib2io._Grib2Message.interpolate` method and 2159 `grib2io.interpolate` function accepts these objects. 2160 """ 2161 gdtn: int 2162 gdt: NDArray 2163 2164 @classmethod 2165 def from_section3(cls, section3): 2166 return cls(section3[4],section3[5:]) 2167 2168 @property 2169 def nx(self): 2170 """Number of grid points in x-direction.""" 2171 return int(self.gdt[7]) 2172 2173 @property 2174 def ny(self): 2175 """Number of grid points in y-direction.""" 2176 return int(self.gdt[8]) 2177 2178 @property 2179 def npoints(self): 2180 """Total number of grid points.""" 2181 return int(self.gdt[7] * self.gdt[8]) 2182 2183 @property 2184 def shape(self): 2185 """Shape of the grid.""" 2186 return (int(self.ny), int(self.nx)) 2187 2188 def to_section3(self): 2189 """Return a full GRIB2 section3 array.""" 2190 return np.array( 2191 [0, self.npoints, 0, 0, self.gdtn] + list(self.gdt) 2192 ).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.
2168 @property 2169 def nx(self): 2170 """Number of grid points in x-direction.""" 2171 return int(self.gdt[7])
Number of grid points in x-direction.
2173 @property 2174 def ny(self): 2175 """Number of grid points in y-direction.""" 2176 return int(self.gdt[8])
Number of grid points in y-direction.
2178 @property 2179 def npoints(self): 2180 """Total number of grid points.""" 2181 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