Using GRIB¶

In this notebook you will see how to:

inspect GRIB data
work with fields and fieldlists
modify fields
convert GRIB to Xarray

Getting the data¶

First we download some GRIB data from the collection of the earthkit example files with from_source. We used the “sample” source for this, but if the data was available as local file on disk we could load it as a “file” source.

import earthkit.data as ekd

d = ekd.from_source("sample", "test4.grib")
d

The returned object provides some basic information but its primary goal is to convert the data into the required representation for further work. The actual data loading is deferred as much as possible, until the data is converted into a given type.

d.available_types

Fieldlists and fields¶

GRIB data can be converted into a FieldList, which represents each GRIB message as a field. In earthkit a field is a horizontal slice of the atmosphere at a given time. In this sense the Field object is generic enough to represent other types of data than GRIB (e.g. NetCDF, GeoTIFF, dictionary data etc.)

# fl is a fieldlist
fl = d.to_fieldlist()

# fl contains 4 fields
len(fl)

ls() lists the fields in the fieldlist.

fl.ls()

# we can iterate through the fields
for f in fl: 
    print(f)

# f is the first field in the fieldlist
f = fl[0]
print(f)

Fields metadata¶

A Field is composed of format independent components, each containing a set of metadata keys. The most important components and keys are as follows:

parameter:
- parameter.variable
- parameter.units
- parameter.chem_variable
time:
- time.base_datetime
- time.valid_datetime
- time.step
vertical:
- vertical.level
- vertical.level_type
geography:
- geography.latitudes
- geography.longitudes
- geography.shape
ensemble:
- ensemble.member

The metadata can be accessed via the get() method.

f.get("time.base_datetime")

We can also use the key methods on the field components.

f.time.base_datetime()

The components contain format independent metadata. This is most obvious when looking at the level type.

# the GRIB metadata would be isobaricInhPa
f.vertical.level_type()

Since the field was created from GRIB the raw (ecCodes) GRIB metadata is also available via the metadata() method.

f.metadata("typeOfLevel")

The get() method also works for raw metadata but the keys must be prefixed with "metadata.".

f.get("metadata.typeOfLevel")

Field units¶

Units are represented by an object that provides normalisation and can used in comparisons (to a string or another Units object).

u = f.parameter.units()
u

u == "kelvin", u == "K"

When possible the Units object is based on a pint.Units object that can be directly accessed when needed.

u_p = u.to_pint()
type(u_p)

Field overview¶

To get an overview about the components/keys of a field simply use the automatic display.

# this is the same as f.describe()
f

Field values¶

Use to_numpy to access the field values as a numpy array. By default it uses the shape of the field.

a = f.to_numpy()
a.shape, a

The values property is also available. Please note this is always return a flattened array.

Field geography¶

# the shape is based on the geography (when available)
f.shape, f.geography.shape()

We can use latlons() to get a tuple of the latitudes and longitudes arrays. Please note “geography.latlons” cannot be used as a key.

f.geography.latlons()

The new way of describing grids is to use a grid spec. It is available as a dict from the field.

f.geography.grid_spec()

Modifying fields¶

Fields can be modified by set(), which generates a new field with updated values/components.

vals = f.values + 2.0
f1 = f.set({"vertical.level": 700, "values": vals})
f1.ls()

# compare the metadata in the old and new fields
f.get("vertical.level"), f1.get("vertical.level")

# compare the values in the old and new fields
f.values.max(), f1.values.max()

Fieldlist subsetting¶

Use sel() to select the fields matching the given metadata conditions.

fl.sel({"parameter.variable": "t"}).ls()

Arithmetics¶

Fields and fieldlists can be used in mathematical expression. The resulting fields/fieldlists will store all the values in memory.

We will compute the geopotential height using fieldlist arithmetic.

# select geopotential fields
z = fl.sel({"parameter.variable": "z"})

# compute the geopotential height
# h and z are fieldlists
h = z / 9.81

# check the results
for zv, hv in zip(z, h):
    print(f"z-max: {zv.values.max()} h-max: {hv.values.max()}")

The metadata in the new “h” fields are still the same as in the “z” fields. We can correct it manually.

h.ls()

h = h.set({"parameter.variable": "gh", "parameter.units": "gpm"})
h.ls()

We repeat the same computation with a field.

# compute the geopotential height
# h and z_first are fields
z_first = z[0]
h = z[0] / 9.81

f"z-max: {z_first.values.max()} h-max: {h.values.max()}"

Converting to Xarray¶

For this exercise we fetch another GRIB file containing more variation in time and level. We check the the available distinct metadata values with the unique() method.

# get GRIB data containing variations in date/time/step and level
d = ekd.from_source("sample", "pl.grib")
fl = d.to_fieldlist()
fl.unique("parameter.variable", "time.base_datetime", "time.step", "vertical.level")

We can convert fieldlists to Xarray using earthkit-data’s own Xarray engine. It comes with a very large number of options, but for our data the defaults are perfectly fine.

ds = fl.to_xarray()
ds

The Xarray we created uses the “forecasts” time dimension mode, with forecast_reference_time and step chosen as the temporal dimensions. The input data also forms a hypercube if we choose the valid datetime as the temporal dimension, so we can use the time_dims option as follows:

ds = fl.to_xarray(time_dims="valid_time")
ds

Xarrays created with earthkit-data have the earthkit accessor. It is an experimental feature and for each DataArray it stores earthkit specific metadata.

ds["t"].earthkit.grid_spec