Earthkit-meteo: interpolate GRIB from model to height levels - ECMWF internal earthkit training, March 2026

In this notebook you will see how to:

interpolate GRIB data from ECMWF model levels to height levels

At the moment, this type of vertical interpolation is only available for array data so we will carry out the computations with the array level interface of earthkit-meteo.

Getting the data¶

The input data contains t (temperature), q (specific humidity), u, v and w fields on model levels 60-137 for one time step. The lnsp (logarithm of surface pressure) and zs (surface geopotential) fields are also available.

import numpy as np
import earthkit.data as ekd
import earthkit.plots as ekp
import earthkit.meteo.vertical.array as vertical

fl = ekd.from_source("sample", "xs_ml_wind_orog.grib").to_fieldlist()
len(fl)

fl[:2].ls()

Preparing the data¶

First, we select the lnsp, t, q and zs fields and extract the values from them as numpy arrays. The t and q fields must be sorted in ascending order with regards to the model levels for the computations.

lnsp = fl.sel({"parameter.variable": "lnsp"})[0]
t = fl.sel({"parameter.variable": "t"}).order_by({"vertical.level": "ascending"})
q = fl.sel({"parameter.variable": "q"}).order_by({"vertical.level": "ascending"})
zs = fl.sel({"parameter.variable": "z", "vertical.level": 1})[0]

sp_vals = np.exp(lnsp.values) # surface pressure
t_vals = t.values
q_vals = q.values
zs_vals = zs .values

Next, we get the hybrid level definitions using earthkit.meteo.vertical.hybrid_level_parameters.

A, B = vertical.hybrid_level_parameters(137)

Computations¶

We carry out the computations in 2 steps:

first, we compute the height on model levels with earthkit.meteo.vertical.height_on_hybrid_level
next, we interpolate the fields to the target heights with earthkit.meteo.vertical.interpolate_monotonic

# compute geometric height above the ground on model levels using numpy arrays
h  = vertical.height_on_hybrid_levels(t_vals, q_vals, zs_vals, A, B, sp_vals, 
                             h_type="geometric", h_reference="ground")

# interpolate t to geometric height (m) levels above the ground
target_h=[1000., 5000.] 

res = vertical.interpolate_monotonic(t_vals, h, target_h, interpolation="linear")

# the result is numpy array
res.shape

Writing the results into a fieldlist¶

Finally, we create a new fieldlist from the results and set the medatadata correctly in it.

t_res = []
for hv, rv in zip(target_h, res):
    f = t[0].set({"vertical.level": hv, "vertical.level_type": "height_above_ground_level", "values": rv})
    t_res.append(f)
t_res = ekd.create_fieldlist(t_res)
t_res.ls()

# check the results visually
ekp.quickplot(t_res)

Interpolating other fields¶

In other to interpolate other vertical fields from the input data we only need to call interpolate_monotonic() since we can reuse the “h” fields. The code below uses a helper method and interpolate the u and v wind components.

def _comp(param):  
    fl_param = fl.sel({"parameter.variable": param}).order_by({"vertical.level": "ascending"})
    vals = fl_param.values
    
    res = vertical.interpolate_monotonic(vals, h, target_h, interpolation="linear")
    
    param_res = []
    for hv, rv in zip(target_h, res):
        f = fl_param[0].set({"vertical.level": hv, "vertical.level_type": "height_above_ground_level", "values": rv})
        param_res.append(f)
    param_res = ekd.create_fieldlist(param_res)
    return param_res

u_res = _comp("u")
v_res = _comp("v")

# check the results visually 
ekp.quickplot(u_res)