Example 2: Reading a GRIB data set¶
GRIB (GRIdded Binary or General Regularly-distributed Information in Binary form) can be read with any of the following libraries or ways:
pygrib
xarray (with the PyNIO engine option)
PyNIO
wgrib2 or CDO to convert to NetCDF
In this example we'll use pygrib to read a WAVEWATCH III dataset, plot the parameter Significant Wave Height (SWH), and select a subregion.
The pygrib library can be installed in several ways: pip, conda (or whatever local package manager you're using for your python installation), or distutils (python setup.py).
- As before, start by importing the required libraries
import matplotlib.pyplot as plt # graphics library
import pygrib
import numpy as np # numeric python library
import cartopy.crs as ccrs # cartographic coord reference system
import cartopy.feature as cfeature # features: land, borders, coastlines
Open the file and display the contents.
GRIB files are record-oriented, and pygrib stores the GRIB records in a list-like object.
You can access individual records by number or just iterate over the list.
gribs = pygrib.open('multi_1.glo_30mext.t00z.f000.grib2')
for grib in gribs:
print(grib)
You can select from the list of GRIB records with long name, short name, or any of the other keys in the file.
Note that the select method returns a list, even for a single message.
swh = gribs(shortName = 'swh',forecastTime = 0)[0] # note single '=", not "=="
- There's a lot of information in each GRIB record
for key in swh.keys(): # doing some python string formatting here
print('{:>25} {}'.format(key, eval('swh.'+key)))
Easiest way to get the data from the GRIB record is to use the
object.data()method.This returns the parameter array and the lat/lon coordinate arrays for the model grid.
data, lats, lons = swh.data()
- Set up the plot figure as before and plot the SWH parameter
plt.figure(dpi = 90)
ax = plt.axes(projection = ccrs.Mercator()) # set up the axis projection
plt.pcolormesh(lons, lats, data, cmap = 'gist_ncar',transform = ccrs.PlateCarree())
cbar = plt.colorbar()
cbar.set_label(swh.shortName+' ('+swh.units+')')
ax.add_feature(cfeature.LAND) # fill in the land areas
ax.coastlines() # use the default low-resolution coastline
gl = ax.gridlines(draw_labels = True) # default is to label all axes.
gl.xlabels_top = False # turn off two of them.
gl.ylabels_right = False
plt.title(swh.name+'\n'+'Valid Date: '+swh.validDate.strftime('%Y%m%d'));
- Now let's look closer at the data itself.
- First off, what is this thing? You can identify the class of any python object with the
type()function.
type(data)
This is a numpy masked array. A masked array is the combination of a standard numpy.ndarray and a Boolean mask of the same shape. The numpy package ensures that the masked data points are excluded from computations.
As before, let's zoom in on Hawaii. We'll use the
object.data()method for that, too. Remember that the "%" is the python symbol for modulo division.
data, lats, lons = swh.data(lon1 = -166%360, lon2 = -150%360,
lat1 = 14, lat2 = 27)
- And now we can proceed to plot the swh parameter.
plt.figure(dpi = 90)
ax = plt.axes(projection = ccrs.Mercator())
plt.contourf(lons, lats, data, transform = ccrs.PlateCarree(), cmap='jet',
levels = np.arange(np.floor(data.min()),np.ceil(data.max()),0.2))
plt.colorbar(label = swh.units)
ax.add_feature(cfeature.LAND,zorder = 1) # zorder will put the land on top
ax.add_feature(cfeature.GSHHSFeature())
gl = ax.gridlines(draw_labels = True)
gl.xlabels_top = False
gl.ylabels_right = False
gl.xlabel_style = {'size': 'x-small'}
gl.ylabel_style = {'size': 'x-small'}
plt.title(swh.name+'\n'+'Valid Date: '+swh.validDate.strftime('%Y%m%d'));
