Loading the data

We load two GRIB files from the earthkit-data sample archive:

• fc_storm_st_jude.grib — the high-resolution deterministic forecast (HRES)

• ens_storm_st_jude.grib — the ensemble forecast (ENS), which contains 51 members: one control forecast (cf, member 0) and 50 perturbed forecasts (pf, members 1–50)

import earthkit.data as ekd

import earthkit.plots as ekp

ds_fc = ekd.from_source("sample", "fc_storm_st_jude.grib")  # hi-res forecast
ds_en = ekd.from_source("sample", "ens_storm_st_jude.grib")  # ensemble forecast

Spaghetti plots

A spaghetti plot overlays a single contour line — drawn at the same value — for every member of an ensemble forecast. The resulting tangle of lines gives an immediate visual impression of forecast uncertainty: where the lines cluster together the ensemble agrees; where they spread apart the forecast is uncertain.

This example uses ECMWF ensemble data from the St Jude storm (October 2013) to produce a spaghetti plot of the 850 hPa geopotential height field.

ds_fc.to_fieldlist().ls()

	parameter.variable	time.valid_datetime	time.base_datetime	time.step	vertical.level	vertical.level_type	ensemble.member	geography.grid_type
0	10fg3	2013-10-28 00:00:00	2013-10-25	3 days 00:00:00	0	surface	0	regular_ll
1	10fg3	2013-10-28 06:00:00	2013-10-25	3 days 06:00:00	0	surface	0	regular_ll
2	10fg3	2013-10-28 12:00:00	2013-10-25	3 days 12:00:00	0	surface	0	regular_ll
3	z	2013-10-28 00:00:00	2013-10-25	3 days 00:00:00	850	pressure	0	regular_ll
4	z	2013-10-28 06:00:00	2013-10-25	3 days 06:00:00	850	pressure	0	regular_ll
5	z	2013-10-28 12:00:00	2013-10-25	3 days 12:00:00	850	pressure	0	regular_ll

Inspecting the data

The deterministic forecast contains several fields. We can see it includes 10-metre wind gust (10fg3) and geopotential height (z) at 850 hPa, each at three lead times (00, 06, and 12 UTC on 28 October 2013).

z_fc = ds_fc.to_fieldlist().sel({"parameter.variable": "z", "vertical.level": 850, "time.step": 78})
z_en = ds_en.to_fieldlist().sel({"parameter.variable": "z", "vertical.level": 850, "time.step": 78})

Selecting a single valid time

We select the 850 hPa geopotential field (z) at step 78 hours from both datasets. This corresponds to 00 UTC on 28 October 2013 — the time when the St Jude storm made landfall over north-west Europe.

Selecting by time.step=78 ensures both the deterministic (z_fc) and ensemble (z_en) fields are at exactly the same valid time, making the comparison meaningful.

ds = z_en.to_xarray()
ds

<xarray.Dataset> Size: 608kB
Dimensions:    (member: 51, latitude: 31, longitude: 48)
Coordinates:
  * member     (member) <U2 408B '0' '1' '10' '11' '12' ... '50' '6' '7' '8' '9'
  * latitude   (latitude) float64 248B 66.0 65.25 64.5 63.75 ... 45.0 44.25 43.5
  * longitude  (longitude) float64 384B -19.5 -18.75 -18.0 ... 14.25 15.0 15.75
Data variables:
    z          (member, latitude, longitude) float64 607kB ...
Attributes:
    Conventions:  CF-1.8
    institution:  ECMWF

xarray.Dataset

• Dimensions:
  • member: 51
  • latitude: 31
  • longitude: 48
• Coordinates: (3)
  • member
    (member)
    <U2
    '0' '1' '10' '11' ... '7' '8' '9'

    standard_name :

    realization

    long_name :

    ensemble member id

    array(['0', '1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19',
           '2', '20', '21', '22', '23', '24', '25', '26', '27', '28', '29', '3',
           '30', '31', '32', '33', '34', '35', '36', '37', '38', '39', '4', '40',
           '41', '42', '43', '44', '45', '46', '47', '48', '49', '5', '50', '6',
           '7', '8', '9'], dtype='<U2')

  • latitude
    (latitude)
    float64
    66.0 65.25 64.5 ... 45.0 44.25 43.5

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    array([66.  , 65.25, 64.5 , 63.75, 63.  , 62.25, 61.5 , 60.75, 60.  , 59.25,
           58.5 , 57.75, 57.  , 56.25, 55.5 , 54.75, 54.  , 53.25, 52.5 , 51.75,
           51.  , 50.25, 49.5 , 48.75, 48.  , 47.25, 46.5 , 45.75, 45.  , 44.25,
           43.5 ])

  • longitude
    (longitude)
    float64
    -19.5 -18.75 -18.0 ... 15.0 15.75

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([-19.5 , -18.75, -18.  , -17.25, -16.5 , -15.75, -15.  , -14.25, -13.5 ,
           -12.75, -12.  , -11.25, -10.5 ,  -9.75,  -9.  ,  -8.25,  -7.5 ,  -6.75,
            -6.  ,  -5.25,  -4.5 ,  -3.75,  -3.  ,  -2.25,  -1.5 ,  -0.75,   0.  ,
             0.75,   1.5 ,   2.25,   3.  ,   3.75,   4.5 ,   5.25,   6.  ,   6.75,
             7.5 ,   8.25,   9.  ,   9.75,  10.5 ,  11.25,  12.  ,  12.75,  13.5 ,
            14.25,  15.  ,  15.75])

• Data variables: (1)
  • z
    (member, latitude, longitude)
    float64
    ...

    standard_name :

    geopotential

    long_name :

    Geopotential

    units :

    meter ** 2 / second ** 2

    level_type :

    pressure

    _earthkit :

    {"message": {"__bytes_b64__": "R1JJQgAAogEAAGqAYo//gIFkA1INChkAAAFOAAAAAAAVAAAAAAAAAAAAAAAAAAAAHgEKBAswMDAxADMAAAAAATMs0QAAAQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAP8AADAAHwEB0IBMLIAAqewAPYYC7gLuAAAAAAAAAAwAgAJEuX1uAAA3Nzc3"}, "bitsPerValue": 16, "grid_spec": {"area": [66, -19.5, 43.5, 15.75], "grid": [0.75, 0.75]}}

    [75888 values with dtype=float64]

• Attributes: (2)

  Conventions :

  CF-1.8

  institution :

  ECMWF

Examining the ensemble structure

Converting to xarray lets us inspect the shape of the ensemble. The dataset has a member dimension of size 51 (members 0–50), with a spatial grid of 31 × 48 points covering roughly 43–66°N, 20°W–16°E — a region centred on north-west Europe.

chart = ekp.Map(figsize=(7, 7))

# the isoline value (geopotential in m²/s²)
cont_level = [12500]

chart.spaghetti(
    z_en,
    levels=cont_level,
    highlight={"metadata.dataType": "cf"},
    label="Ensemble members",
)

chart.land()
chart.coastlines()
chart.borders()
chart.gridlines()

chart.title(
    "ECMWF Run: {base_time:%Y-%m-%d %H} UTC (+{lead_time}h) {variable_name} {level} hPa",
)

chart.legend()

chart.show()

Plotting the spaghetti

chart.spaghetti() iterates over every member in the ensemble and draws a contour line at the chosen level for each one.

Key parameters:

• ``levels=[12500]`` — the geopotential value (in m²/s²) at which to draw the contour. A single isoline is drawn per member, which is the defining characteristic of a spaghetti plot.

• ``highlight={“metadata.dataType”: “cf”}`` — picks out the control forecast member (member 0) and renders it in a distinct style so it stands apart from the perturbed members.

• ``label=”foo”`` — assigns a legend label to the ensemble member lines.

Where the 51 lines overlap closely the ensemble is confident about the position of that height contour. Where they diverge — most visibly around the storm centre — the forecast is uncertain.