Vertical (Xarray): interpolating from pressure to pressure levels¶

This notebook demonstrates how to interpolate Xarray data on pressure levels to another set of pressure levels with interpolate_monotonic(). We will use earthkit-data for getting the input data and earthkit-plots for plotting.

[1]:

import earthkit.data as ekd
import earthkit.plots as ekp

import earthkit.meteo.vertical.xarray as vertical

Getting the data¶

The input is NetCDF data on pressure levels for multiple times. We fetch the data and load it into an Xarray.

[2]:

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

[2]:

<xarray.Dataset> Size: 35MB
Dimensions:         (valid_time: 4, pressure_level: 3, latitude: 301,
                     longitude: 601)
Coordinates:
  * valid_time      (valid_time) datetime64[ns] 32B 2016-09-26 ... 2016-09-26...
  * pressure_level  (pressure_level) float64 24B 925.0 850.0 700.0
  * latitude        (latitude) float64 2kB 90.0 89.75 89.5 ... 15.5 15.25 15.0
  * longitude       (longitude) float64 5kB -100.0 -99.75 -99.5 ... 49.75 50.0
Data variables:
    q               (valid_time, pressure_level, latitude, longitude) float32 9MB ...
    t               (valid_time, pressure_level, latitude, longitude) float32 9MB ...
    u               (valid_time, pressure_level, latitude, longitude) float32 9MB ...
    v               (valid_time, pressure_level, latitude, longitude) float32 9MB ...
Attributes:
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2026-03-12T10:30 GRIB to CDM+CF via cfgrib-0.9.1...

xarray.Dataset

Dimensions:
- valid_time: 4
- pressure_level: 3
- latitude: 301
- longitude: 601
Coordinates: (4)
- valid_time
  (valid_time)
  datetime64[ns]
  2016-09-26 ... 2016-09-26T18:00:00
  long_name :
  time
  standard_name :
  time
```
array(['2016-09-26T00:00:00.000000000', '2016-09-26T06:00:00.000000000',
       '2016-09-26T12:00:00.000000000', '2016-09-26T18:00:00.000000000'],
      dtype='datetime64[ns]')
```
- pressure_level
  (pressure_level)
  float64
  925.0 850.0 700.0
  long_name :
  pressure
  units :
  hPa
  positive :
  down
  stored_direction :
  decreasing
  standard_name :
  air_pressure
```
array([925., 850., 700.])
```
- latitude
  (latitude)
  float64
  90.0 89.75 89.5 ... 15.5 15.25 15.0
  units :
  degrees_north
  standard_name :
  latitude
  long_name :
  latitude
  stored_direction :
  decreasing
```
array([90.  , 89.75, 89.5 , ..., 15.5 , 15.25, 15.  ], shape=(301,))
```
- longitude
  (longitude)
  float64
  -100.0 -99.75 -99.5 ... 49.75 50.0
  units :
  degrees_east
  standard_name :
  longitude
  long_name :
  longitude
```
array([-100.  ,  -99.75,  -99.5 , ...,   49.5 ,   49.75,   50.  ], shape=(601,))
```
Data variables: (4)
- q
  (valid_time, pressure_level, latitude, longitude)
  float32
  ...
  GRIB_paramId :
  133
  GRIB_dataType :
  an
  GRIB_numberOfPoints :
  180901
  GRIB_typeOfLevel :
  isobaricInhPa
  GRIB_stepUnits :
  1
  GRIB_stepType :
  instant
  GRIB_gridType :
  regular_ll
  GRIB_uvRelativeToGrid :
  0
  GRIB_NV :
  0
  GRIB_Nx :
  601
  GRIB_Ny :
  301
  GRIB_cfName :
  specific_humidity
  GRIB_cfVarName :
  q
  GRIB_gridDefinitionDescription :
  Latitude/Longitude Grid
  GRIB_iDirectionIncrementInDegrees :
  0.25
  GRIB_iScansNegatively :
  0
  GRIB_jDirectionIncrementInDegrees :
  0.25
  GRIB_jPointsAreConsecutive :
  0
  GRIB_jScansPositively :
  0
  GRIB_latitudeOfFirstGridPointInDegrees :
  90.0
  GRIB_latitudeOfLastGridPointInDegrees :
  15.0
  GRIB_longitudeOfFirstGridPointInDegrees :
  -100.0
  GRIB_longitudeOfLastGridPointInDegrees :
  50.0
  GRIB_missingValue :
  3.4028234663852886e+38
  GRIB_name :
  Specific humidity
  GRIB_shortName :
  q
  GRIB_totalNumber :
  0
  GRIB_units :
  kg kg**-1
  long_name :
  Specific humidity
  units :
  kg kg**-1
  standard_name :
  specific_humidity
```
[2170812 values with dtype=float32]
```
- t
  (valid_time, pressure_level, latitude, longitude)
  float32
  ...
  GRIB_paramId :
  130
  GRIB_dataType :
  an
  GRIB_numberOfPoints :
  180901
  GRIB_typeOfLevel :
  isobaricInhPa
  GRIB_stepUnits :
  1
  GRIB_stepType :
  instant
  GRIB_gridType :
  regular_ll
  GRIB_uvRelativeToGrid :
  0
  GRIB_NV :
  0
  GRIB_Nx :
  601
  GRIB_Ny :
  301
  GRIB_cfName :
  air_temperature
  GRIB_cfVarName :
  t
  GRIB_gridDefinitionDescription :
  Latitude/Longitude Grid
  GRIB_iDirectionIncrementInDegrees :
  0.25
  GRIB_iScansNegatively :
  0
  GRIB_jDirectionIncrementInDegrees :
  0.25
  GRIB_jPointsAreConsecutive :
  0
  GRIB_jScansPositively :
  0
  GRIB_latitudeOfFirstGridPointInDegrees :
  90.0
  GRIB_latitudeOfLastGridPointInDegrees :
  15.0
  GRIB_longitudeOfFirstGridPointInDegrees :
  -100.0
  GRIB_longitudeOfLastGridPointInDegrees :
  50.0
  GRIB_missingValue :
  3.4028234663852886e+38
  GRIB_name :
  Temperature
  GRIB_shortName :
  t
  GRIB_totalNumber :
  0
  GRIB_units :
  K
  long_name :
  Temperature
  units :
  K
  standard_name :
  air_temperature
```
[2170812 values with dtype=float32]
```
- u
  (valid_time, pressure_level, latitude, longitude)
  float32
  ...
  GRIB_paramId :
  131
  GRIB_dataType :
  an
  GRIB_numberOfPoints :
  180901
  GRIB_typeOfLevel :
  isobaricInhPa
  GRIB_stepUnits :
  1
  GRIB_stepType :
  instant
  GRIB_gridType :
  regular_ll
  GRIB_uvRelativeToGrid :
  0
  GRIB_NV :
  0
  GRIB_Nx :
  601
  GRIB_Ny :
  301
  GRIB_cfName :
  eastward_wind
  GRIB_cfVarName :
  u
  GRIB_gridDefinitionDescription :
  Latitude/Longitude Grid
  GRIB_iDirectionIncrementInDegrees :
  0.25
  GRIB_iScansNegatively :
  0
  GRIB_jDirectionIncrementInDegrees :
  0.25
  GRIB_jPointsAreConsecutive :
  0
  GRIB_jScansPositively :
  0
  GRIB_latitudeOfFirstGridPointInDegrees :
  90.0
  GRIB_latitudeOfLastGridPointInDegrees :
  15.0
  GRIB_longitudeOfFirstGridPointInDegrees :
  -100.0
  GRIB_longitudeOfLastGridPointInDegrees :
  50.0
  GRIB_missingValue :
  3.4028234663852886e+38
  GRIB_name :
  U component of wind
  GRIB_shortName :
  u
  GRIB_totalNumber :
  0
  GRIB_units :
  m s**-1
  long_name :
  U component of wind
  units :
  m s**-1
  standard_name :
  eastward_wind
```
[2170812 values with dtype=float32]
```
- v
  (valid_time, pressure_level, latitude, longitude)
  float32
  ...
  GRIB_paramId :
  132
  GRIB_dataType :
  an
  GRIB_numberOfPoints :
  180901
  GRIB_typeOfLevel :
  isobaricInhPa
  GRIB_stepUnits :
  1
  GRIB_stepType :
  instant
  GRIB_gridType :
  regular_ll
  GRIB_uvRelativeToGrid :
  0
  GRIB_NV :
  0
  GRIB_Nx :
  601
  GRIB_Ny :
  301
  GRIB_cfName :
  northward_wind
  GRIB_cfVarName :
  v
  GRIB_gridDefinitionDescription :
  Latitude/Longitude Grid
  GRIB_iDirectionIncrementInDegrees :
  0.25
  GRIB_iScansNegatively :
  0
  GRIB_jDirectionIncrementInDegrees :
  0.25
  GRIB_jPointsAreConsecutive :
  0
  GRIB_jScansPositively :
  0
  GRIB_latitudeOfFirstGridPointInDegrees :
  90.0
  GRIB_latitudeOfLastGridPointInDegrees :
  15.0
  GRIB_longitudeOfFirstGridPointInDegrees :
  -100.0
  GRIB_longitudeOfLastGridPointInDegrees :
  50.0
  GRIB_missingValue :
  3.4028234663852886e+38
  GRIB_name :
  V component of wind
  GRIB_shortName :
  v
  GRIB_totalNumber :
  0
  GRIB_units :
  m s**-1
  long_name :
  V component of wind
  units :
  m s**-1
  standard_name :
  northward_wind
```
[2170812 values with dtype=float32]
```
Attributes: (6)
GRIB_centre :
ecmf
GRIB_centreDescription :
European Centre for Medium-Range Weather Forecasts
GRIB_subCentre :
0
Conventions :
CF-1.7
institution :
European Centre for Medium-Range Weather Forecasts
history :
2026-03-12T10:30 GRIB to CDM+CF via cfgrib-0.9.15.1/ecCodes-2.42.0 with {"source": "tmpfsay9t6q/data.grib", "filter_by_keys": {"stream": ["oper"], "stepType": ["instant"]}, "encode_cf": ["parameter", "time", "geography", "vertical"]}

Using interpolate_monotonic¶

We use interpolate_monotonic() to carry out the interpolation.

interpolate_monotonic() uses the coords, target_coords arguments to specify the source and target coordinates.

[3]:

# get dataarray for temperature on pressure levels
t = ds.t

# define the pressure coordinates of t (in Pa)
coords = t.pressure_level * 100  # hPa to Pa

# define the target pressure coordinates
target_coords = [87000.0, 90000]  # Pa

t_res = vertical.interpolate_monotonic(
    t,  # data to interpolate
    coords=t.pressure_level * 100,
    target_coords=target_coords,
    interpolation="log",
    vertical_dim="pressure_level",
)

# dump the result
t_res

[3]:

<xarray.DataArray 't' (pressure_level: 2, valid_time: 4, latitude: 301,
                       longitude: 601)> Size: 6MB
array([[[[268.80438, 268.80438, 268.80438, ..., 268.80438, 268.80438,
          268.80438],
         [268.96716, 268.97284, 268.97708, ..., 268.8706 , 268.86777,
          268.86636],
         [269.39813, 269.4061 , 269.41602, ..., 269.07675, 269.0725 ,
          269.06628],
         ...,
         [292.64767, 292.67758, 292.7104 , ..., 299.59018, 299.24353,
          298.86963],
         [292.65625, 292.6461 , 292.67627, ..., 299.43314, 298.91785,
          298.5941 ],
         [292.65518, 292.59857, 292.64688, ..., 298.97372, 298.7429 ,
          298.5587 ]],

        [[267.9577 , 267.9577 , 267.9577 , ..., 267.9577 , 267.9577 ,
          267.9577 ],
         [267.7815 , 267.78717, 267.7923 , ..., 268.49957, 268.49673,
          268.4925 ],
         [267.50983, 267.52063, 267.53088, ..., 268.8865 , 268.87997,
          268.8743 ],
...
         [292.8984 , 292.83447, 292.9217 , ..., 305.53693, 305.3479 ,
          305.1117 ],
         [292.96344, 292.744  , 292.78485, ..., 303.21362, 303.451  ,
          302.9143 ],
         [293.06302, 292.83704, 292.7677 , ..., 298.51385, 299.11368,
          299.00436]],

        [[267.94824, 267.94824, 267.94824, ..., 267.94824, 267.94824,
          267.94824],
         [267.95895, 267.9596 , 267.96088, ..., 268.09116, 268.0899 ,
          268.08926],
         [268.01514, 268.0164 , 268.01767, ..., 268.24216, 268.2409 ,
          268.23962],
         ...,
         [293.54477, 293.45462, 293.47855, ..., 301.37518, 301.30316,
          300.94415],
         [293.44983, 293.38193, 293.34415, ..., 298.6285 , 298.39713,
          298.66498],
         [293.35663, 293.3362 , 293.32257, ..., 297.88422, 297.91034,
          298.24277]]]], shape=(2, 4, 301, 601), dtype=float32)
Coordinates:
  * pressure_level  (pressure_level) float64 16B 8.7e+04 9e+04
  * valid_time      (valid_time) datetime64[ns] 32B 2016-09-26 ... 2016-09-26...
  * latitude        (latitude) float64 2kB 90.0 89.75 89.5 ... 15.5 15.25 15.0
  * longitude       (longitude) float64 5kB -100.0 -99.75 -99.5 ... 49.75 50.0
Attributes: (12/31)
    GRIB_paramId:                             130
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      180901
    GRIB_typeOfLevel:                         isobaricInhPa
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    ...                                       ...
    GRIB_shortName:                           t
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                Temperature
    units:                                    K
    standard_name:                            air_temperature

xarray.DataArray

't'

pressure_level: 2
valid_time: 4
latitude: 301
longitude: 601

268.8 268.8 268.8 268.8 268.8 268.8 ... 298.3 297.9 297.9 297.9 298.2

array([[[[268.80438, 268.80438, 268.80438, ..., 268.80438, 268.80438,
          268.80438],
         [268.96716, 268.97284, 268.97708, ..., 268.8706 , 268.86777,
          268.86636],
         [269.39813, 269.4061 , 269.41602, ..., 269.07675, 269.0725 ,
          269.06628],
         ...,
         [292.64767, 292.67758, 292.7104 , ..., 299.59018, 299.24353,
          298.86963],
         [292.65625, 292.6461 , 292.67627, ..., 299.43314, 298.91785,
          298.5941 ],
         [292.65518, 292.59857, 292.64688, ..., 298.97372, 298.7429 ,
          298.5587 ]],

        [[267.9577 , 267.9577 , 267.9577 , ..., 267.9577 , 267.9577 ,
          267.9577 ],
         [267.7815 , 267.78717, 267.7923 , ..., 268.49957, 268.49673,
          268.4925 ],
         [267.50983, 267.52063, 267.53088, ..., 268.8865 , 268.87997,
          268.8743 ],
...
         [292.8984 , 292.83447, 292.9217 , ..., 305.53693, 305.3479 ,
          305.1117 ],
         [292.96344, 292.744  , 292.78485, ..., 303.21362, 303.451  ,
          302.9143 ],
         [293.06302, 292.83704, 292.7677 , ..., 298.51385, 299.11368,
          299.00436]],

        [[267.94824, 267.94824, 267.94824, ..., 267.94824, 267.94824,
          267.94824],
         [267.95895, 267.9596 , 267.96088, ..., 268.09116, 268.0899 ,
          268.08926],
         [268.01514, 268.0164 , 268.01767, ..., 268.24216, 268.2409 ,
          268.23962],
         ...,
         [293.54477, 293.45462, 293.47855, ..., 301.37518, 301.30316,
          300.94415],
         [293.44983, 293.38193, 293.34415, ..., 298.6285 , 298.39713,
          298.66498],
         [293.35663, 293.3362 , 293.32257, ..., 297.88422, 297.91034,
          298.24277]]]], shape=(2, 4, 301, 601), dtype=float32)

Coordinates: (4)
- pressure_level
  (pressure_level)
  float64
  8.7e+04 9e+04
```
array([87000., 90000.])
```
- valid_time
  (valid_time)
  datetime64[ns]
  2016-09-26 ... 2016-09-26T18:00:00
  long_name :
  time
  standard_name :
  time
```
array(['2016-09-26T00:00:00.000000000', '2016-09-26T06:00:00.000000000',
       '2016-09-26T12:00:00.000000000', '2016-09-26T18:00:00.000000000'],
      dtype='datetime64[ns]')
```
- latitude
  (latitude)
  float64
  90.0 89.75 89.5 ... 15.5 15.25 15.0
  units :
  degrees_north
  standard_name :
  latitude
  long_name :
  latitude
  stored_direction :
  decreasing
```
array([90.  , 89.75, 89.5 , ..., 15.5 , 15.25, 15.  ], shape=(301,))
```
- longitude
  (longitude)
  float64
  -100.0 -99.75 -99.5 ... 49.75 50.0
  units :
  degrees_east
  standard_name :
  longitude
  long_name :
  longitude
```
array([-100.  ,  -99.75,  -99.5 , ...,   49.5 ,   49.75,   50.  ], shape=(601,))
```
Attributes: (31)
GRIB_paramId :
130
GRIB_dataType :
an
GRIB_numberOfPoints :
180901
GRIB_typeOfLevel :
isobaricInhPa
GRIB_stepUnits :
1
GRIB_stepType :
instant
GRIB_gridType :
regular_ll
GRIB_uvRelativeToGrid :
0
GRIB_NV :
0
GRIB_Nx :
601
GRIB_Ny :
301
GRIB_cfName :
air_temperature
GRIB_cfVarName :
t
GRIB_gridDefinitionDescription :
Latitude/Longitude Grid
GRIB_iDirectionIncrementInDegrees :
0.25
GRIB_iScansNegatively :
0
GRIB_jDirectionIncrementInDegrees :
0.25
GRIB_jPointsAreConsecutive :
0
GRIB_jScansPositively :
0
GRIB_latitudeOfFirstGridPointInDegrees :
90.0
GRIB_latitudeOfLastGridPointInDegrees :
15.0
GRIB_longitudeOfFirstGridPointInDegrees :
-100.0
GRIB_longitudeOfLastGridPointInDegrees :
50.0
GRIB_missingValue :
3.4028234663852886e+38
GRIB_name :
Temperature
GRIB_shortName :
t
GRIB_totalNumber :
0
GRIB_units :
K
long_name :
Temperature
units :
K
standard_name :
air_temperature

Xarray based interpolation is still under development and not all the metadata (e.g. level units as you can see in the dump above) is set on the result.

We plot a slice of the results.

[4]:

ekp.quickplot(t_res.isel(valid_time=0))

[4]:

<earthkit.plots.components.figures.Figure at 0x135669160>

../../_images/tutorials_vertical_interpolate_pl_to_pl_xarray_11_1.png

[ ]: