Interpolating HEALPix GRIB fields

This example shows how to interpolate GRIB data defined on a HEALPix nested grid. We will also see how to inspect and plot the resulting data and how to convert it to xarray.

To make this notebook work earthkit-data (version >= 0.6.0) and earthkit-plots have to be installed. The data will be represented as an earthkit-data GRIB FieldList.

Interpolation

from earthkit.regrid import interpolate
from earthkit.data import from_source

# Get HEALPix nested GRIB data containing two fields.
ds = from_source(
    "url",
    "https://sites.ecmwf.int/repository/earthkit-regrid/examples/H8_nested_multi.grib2")

# the target grid is a global 5x5 degree regular latitude-longitude grid
out_grid = {"grid": [5,5]}

# perform interpolation for each field and add results
# to a new fieldlist stored in memory
r = interpolate(ds, out_grid=out_grid, method="linear")

d = r.data()
lat = d[0]
lon = d[1]
vals = d[2:]
lat.shape, lon.shape, vals.shape

/var/folders/93/w0p869rx17q98wxk83gn9ys40000gn/T/ipykernel_42625/2538999610.py:14: DeprecatedWarning: interpolate is deprecated as of 0.5.0. Use regrid() instead
  r = interpolate(ds, out_grid=out_grid, method="linear")

((37, 72), (37, 72), (2, 37, 72))

Please note that the interpolation only works between a predefined set of global grids. See the Pre-generated matrices for details.

Plotting the results

import earthkit.plots

figure = earthkit.plots.Figure(size=(9, 6), rows=1, columns=2)
figure.add_map().plot(r[0])
figure.add_map().plot(r[1])
figure.coastlines()
figure.subplot_titles()
figure.legend()
figure.gridlines(sharey=True)
figure.show()

/Users/cgr/venv/mir/lib/python3.13/site-packages/earthkit/plots/components/subplots.py:727: UserWarning: `plot` is deprecated. Use `quickplot` instead.
  warnings.warn("`plot` is deprecated. Use `quickplot` instead.")
/Users/cgr/venv/mir/lib/python3.13/site-packages/earthkit/plots/components/subplots.py:727: UserWarning: `plot` is deprecated. Use `quickplot` instead.
  warnings.warn("`plot` is deprecated. Use `quickplot` instead.")

Converting the results to xarray

r.to_xarray()

<xarray.Dataset> Size: 44kB
Dimensions:    (step: 2, latitude: 37, longitude: 72)
Coordinates:
  * step       (step) timedelta64[ns] 16B 00:00:00 12:00:00
  * latitude   (latitude) float64 296B 90.0 85.0 80.0 75.0 ... -80.0 -85.0 -90.0
  * longitude  (longitude) float64 576B 0.0 5.0 10.0 15.0 ... 345.0 350.0 355.0
Data variables:
    2t         (step, latitude, longitude) float64 43kB ...
Attributes:
    param:        2t
    paramId:      167
    class:        od
    stream:       oper
    levtype:      sfc
    type:         fc
    expver:       0001
    date:         20240323
    time:         1200
    domain:       g
    Conventions:  CF-1.8
    institution:  ECMWF

xarray.Dataset

• Dimensions:
  • step: 2
  • latitude: 37
  • longitude: 72
• Coordinates: (3)
  • step
    (step)
    timedelta64[ns]
    00:00:00 12:00:00

    array([             0, 43200000000000], dtype='timedelta64[ns]')

  • latitude
    (latitude)
    float64
    90.0 85.0 80.0 ... -85.0 -90.0

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    array([ 90.,  85.,  80.,  75.,  70.,  65.,  60.,  55.,  50.,  45.,  40.,  35.,
            30.,  25.,  20.,  15.,  10.,   5.,   0.,  -5., -10., -15., -20., -25.,
           -30., -35., -40., -45., -50., -55., -60., -65., -70., -75., -80., -85.,
           -90.])

  • longitude
    (longitude)
    float64
    0.0 5.0 10.0 ... 345.0 350.0 355.0

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([  0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,  55.,
            60.,  65.,  70.,  75.,  80.,  85.,  90.,  95., 100., 105., 110., 115.,
           120., 125., 130., 135., 140., 145., 150., 155., 160., 165., 170., 175.,
           180., 185., 190., 195., 200., 205., 210., 215., 220., 225., 230., 235.,
           240., 245., 250., 255., 260., 265., 270., 275., 280., 285., 290., 295.,
           300., 305., 310., 315., 320., 325., 330., 335., 340., 345., 350., 355.])

• Data variables: (1)
  • 2t
    (step, latitude, longitude)
    float64
    ...

    standard_name :

    air_temperature

    long_name :

    2 metre temperature

    units :

    K

    _earthkit :

    {'message': b'GRIB\xff\xff\x00\x02\x00\x00\x00\x00\x00\x00\x00\xc4\x00\x00\x00\x15\x01\x00b\x00\x00 \x00\x01\x07\xe8\x03\x17\x0c\x00\x00\x00\x01\x00\x00\x00\x11\x02\x00\x01\x00\x01\x00\t\x04\x010001\x00\x00\x00H\x03\x00\x00\x00\nh\x00\x00\x00\x00\x00\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\x00\x00\x00H\x00\x00\x00%\x00\x00\x00\x00\xff\xff\xff\xff\x05]J\x80\x00\x00\x00\x000\x85]J\x80\x15(\xde\xc0\x00LK@\x00LK@\x00\x00\x00\x00"\x04\x00\x00\x00\x00\x00\x00\x00\xff\x9a\xff\xff\xff\x01\x00\x00\x00\x00g\x00\x00\x00\x00\x02\xff\xff\xff\xff\xff\xff\x00\x00\x00\x15\x05\x00\x00\nh\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x06\x06\xff\x00\x00\x00\x05\x077777', 'bitsPerValue': 24}

    [5328 values with dtype=float64]

• Indexes: (3)
  • step
    PandasIndex

    PandasIndex(TimedeltaIndex(['0 days 00:00:00', '0 days 12:00:00'], dtype='timedelta64[ns]', name='step', freq=None))

  • latitude
    PandasIndex

    PandasIndex(Index([ 90.0,  85.0,  80.0,  75.0,  70.0,  65.0,  60.0,  55.0,  50.0,  45.0,
            40.0,  35.0,  30.0,  25.0,  20.0,  15.0,  10.0,   5.0,   0.0,  -5.0,
           -10.0, -15.0, -20.0, -25.0, -30.0, -35.0, -40.0, -45.0, -50.0, -55.0,
           -60.0, -65.0, -70.0, -75.0, -80.0, -85.0, -90.0],
          dtype='float64', name='latitude'))

  • longitude
    PandasIndex

    PandasIndex(Index([  0.0,   5.0,  10.0,  15.0,  20.0,  25.0,  30.0,  35.0,  40.0,  45.0,
            50.0,  55.0,  60.0,  65.0,  70.0,  75.0,  80.0,  85.0,  90.0,  95.0,
           100.0, 105.0, 110.0, 115.0, 120.0, 125.0, 130.0, 135.0, 140.0, 145.0,
           150.0, 155.0, 160.0, 165.0, 170.0, 175.0, 180.0, 185.0, 190.0, 195.0,
           200.0, 205.0, 210.0, 215.0, 220.0, 225.0, 230.0, 235.0, 240.0, 245.0,
           250.0, 255.0, 260.0, 265.0, 270.0, 275.0, 280.0, 285.0, 290.0, 295.0,
           300.0, 305.0, 310.0, 315.0, 320.0, 325.0, 330.0, 335.0, 340.0, 345.0,
           350.0, 355.0],
          dtype='float64', name='longitude'))

• Attributes: (12)

  param :

  2t

  paramId :

  167

  class :

  od

  stream :

  oper

  levtype :

  sfc

  type :

  fc

  expver :

  0001

  date :

  20240323

  time :

  1200

  domain :

  g

  Conventions :

  CF-1.8

  institution :

  ECMWF

Writing the results to disk

Write the resulting fieldlist to disk:

out_file = "_res_H8_nested_to_5x5.grib"
r.save(out_file)

/var/folders/93/w0p869rx17q98wxk83gn9ys40000gn/T/ipykernel_42625/2197547179.py:2: DeprecatedWarning: save is deprecated as of 0.13.0. Use to_target() instead
  r.save(out_file)