Earthkit-meteo: thermodynamics - ECMWF internal earthkit training, March 2026

In this notebook you will see how to:

compute the potential temperature from NetCDF/Xarray data
compute the relative humidity from NetCDF/Xarray data

Getting the input data¶

First, we get some NetCDF data.

import earthkit.data as ekd
from earthkit.meteo import thermo
import earthkit.plots as ekp

ds = ekd.from_source("sample", "era5_tquv_pl_subarea.nc").to_xarray()
ds

Potential temperature¶

We compute the potential temperature with earthkit.meteo.thermo.potential_temperature. The pressure is required in Pa units, but it is not automatically scaled for us. Since the pressure in the data is in hPa we have to scale it manually.

# input: temperature (K) and pressure (Pa)
theta = thermo.potential_temperature(ds.t, ds.pressure_level*100.)
theta

Earthkit-meteo is still work in progress and the metadata is not yet set in the resulting DataArray. We will set it manually.

theta.attrs["standard_name"] = "air_potential_temperature"
theta.attrs["long_name"] = "Potential temperature"
theta = theta.rename("pt")

# select one level and step and plot input and resulting fields
t1 = ds.t.sel(pressure_level=850, valid_time="2016-09-26T12")
theta1 = theta.sel(pressure_level=850, valid_time="2016-09-26T12")

fig = ekp.Figure(rows=1, columns=2)

fig.add_map().contourf(t1, cmap="plasma")
fig.add_map().contourf(theta1, cmap="viridis")

fig.coastlines()
fig.subplot_titles()
fig.legend()

fig.show()

Relative humidity¶

We compute the relative humidity with earthkit.meteo.thermo.relative_humidity_from_specific_humidity. The pressure requires the same treatment as above.

# input: temperature (K), specific humidity (kg/kg) and pressure (Pa)
r = thermo.relative_humidity_from_specific_humidity(ds.t,ds.q,ds.pressure_level*100.)
r

We set the resulting metadata.

r.attrs["standard_name"] = "relative_humidity"
r.attrs["long_name"] = "Relative humidity"
r = r.rename("r")

# plot one level and step from the results
r1 = r.sel(pressure_level=850, valid_time="2016-09-26T12")
ekp.quickplot(r1)