Heavy Rain - AIFS - Open Data Visio

In the beginning of February 2025, Northern Queensland was hit by extreme rainfall event and large areas of land were flooded. Australia’s Bureau of Meteorology reported that the area had received six months of rainfall in three days.

The example shows analysis of the following parameters:

msl mean sea level pressure, and
tp total precipitation of the AIFS datasets on 31 January at 00 UTC in North Queensland (20.77° S, 144.79° E), Australia.

1. Set Up Your Environment and Find ECMWF Open Data¶

Open data will be downloaded from a publicly available Amazon S3 Bucket. First, the following Python libraries need to be installed in the current Jupyter kernel:

ecmwf-opendata to download data and
earthkit to analyse and plot the data.

If the packages are not installed yet, uncomment the code below and run it.

# !pip3 install earthkit ecmwf-opendata

from ecmwf.opendata import Client
import earthkit.data as ekd
import earthkit.plots as ekp

List of parameters to retrieve from open datasets¶

The selected values below can be modified.

Parameters available on a single level:

PARAM_SFC = ["msl", "tp"]
LEVELTYPE = "sfc"

DATES = [20250130, 20250131]
TIME = 0
STEPS = 12
STREAM = "oper"
TYPE = "fc"
MODEL = "aifs"

Get the data using the ECMWF Open Data API¶

def get_open_data(date, time, step, stream, _type, model, param, leveltype, levelist=[]):
    client = Client(source="aws")
    list_of_files = []
    # Get the data for all dates
    for _date in DATES:
        filename = f"{model}_{''.join(param)}_{''.join(map(str, levelist))}_{_date}.grib2" if levelist else f"{model}_{''.join(param)}_{leveltype}_{_date}.grib2"
        data = client.retrieve(
            date=_date,
            time=time,
            step=step,
            stream=stream,
            type=_type,
            levtype=leveltype,
            levelist=levelist,
            param=param,
            model=model,
            target=filename
            )
        list_of_files.append(filename)
    return data, list_of_files

2. Total precipitation and mean sea level pressure¶

data, list_of_files = get_open_data(date=DATES,
                                    time=TIME,
                                    step=STEPS,
                                    stream=STREAM,
                                    _type=TYPE,
                                    model=MODEL,
                                    param=PARAM_SFC,
                                    leveltype=LEVELTYPE,
                                    levelist=[])

# Select AIFS model data on 31 January 2025
ds = ekd.from_source("file", list_of_files[1])
ds.ls()

When using the head() method, a selected number of rows n and information about the fields in the file we downloaded will be displayed.

ds.head(n=2, namespace="parameter")

A unit of total precipitation is kg/m $^2$ . 1 kg of rainwater fills an area of 1 m $^2$ with the water of height 1 mm.
In the ECMWF Open Charts, total precipitation is also plotted in millimetres.

tp = ds.sel(shortName="tp")

We will plot mean sea level pressure data in hPa, therefore we need to divide them by 100.

msl = ds.sel(shortName="msl")
md_msl = msl.metadata()

msl_hPa = msl.values / 100
msl_3101 = ekd.FieldList.from_array(msl_hPa, md_msl)
msl_3101.ls()

3. Total precipitation¶

The input values can be set here.

PARAM_SFC = "tp"
LEVELS = []
LEVELTYPE = "sfc"

DATES = 20250130
TIME = 0
STEPS = [12, 24, 36, 48, 60, 72, 84, 96]
STREAM = "oper"
TYPE = "fc"
MODEL = "aifs"

The function below is similar to the one defined above.

filename = f"{MODEL}_{PARAM_SFC}_{DATES}_{'-'.join(map(str, STEPS))}.grib2"

def get_open_data_Nsteps(filename, date, time, step, stream, _type, model, param, leveltype, levelist=[]):
    client = Client(source="aws")
    # Get the data for all steps
    data = client.retrieve(
        date=date,
        time=time,
        step=step,
        stream=stream,
        type=_type,
        levtype=leveltype,
        levelist=levelist,
        param=param,
        model=model,
        target=filename
        )
    return data

data = get_open_data_Nsteps(filename,
                            date=DATES,
                            time=TIME,
                            step=STEPS,
                            stream=STREAM,
                            _type=TYPE,
                            model=MODEL,
                            param=PARAM_SFC,
                            leveltype=LEVELTYPE,
                            levelist=[])

ds_tp = ekd.from_source("file", filename)
ds_tp.ls()

The tp parameter gives information about total accumulated rainfall from the start of the forecast onwards. For instance, step=12 indicates accumulated precipitation from 00 UTC until 12 UTC, step=96 from 00 UTC up to 4 days ahead.

4. Data visualisation¶

The plot below shows analysis of mean sea level pressure and total precipitation on 31 January 2025.

chart = ekp.Map(domain="Australia")

tp_style = ekp.styles.Style(
    colors="Spectral_r",
    levels=range(1, 250, 5),
    units="mm",
)

chart.grid_cells(tp, style=tp_style)
chart.contour(msl_3101, legend_style="None")

chart.coastlines(resolution="low")
chart.gridlines()
chart.cities(adjust_labels=True)

chart.legend(location="bottom", label="{variable_name} ({units})")

chart.title(
    "AIFS: {variable_name} over {domain}\n"
    "{base_time:%Y-%m-%d %H} UTC (+{lead_time}h)\n",
    fontsize=14, horizontalalignment="center",
)
chart.save(f"./plots/{''.join(PARAM_SFC)}_{MODEL}_{DATES[1]}{TIME}-{STEPS}h.png")
chart.show()

The plots below show analyses of total precipitation from 30 January at 00 UTC to 3 February, every 12 hours.

figure = ekp.Figure(domain=[130, 155, -35, -10], size=(9, 8), rows=3, columns=4)

dstp_style = ekp.styles.Style(
    colors="Spectral_r",
    levels=range(1, 500, 5),
    units="mm",
)

for i in range(8):
    figure.add_map(1+i//4, i%4)

figure.contourf(ds_tp, style=dstp_style)

figure.coastlines()
figure.gridlines()

figure.legend(label="{variable_name} ({units})")

figure.subplot_titles("{base_time:%Y-%m-%d %H} UTC (+{lead_time}h)")
figure.title(
    "AIFS: {base_time:%Y-%m-%d %H} UTC\n {variable_name} over Northern Queensland\n\n",
    fontsize=14, horizontalalignment="center",
)
figure.save(fname=f"./plots/{PARAM_SFC}_{MODEL}_{DATES}{TIME}-{'-'.join(map(str, STEPS))}h.png")
figure.show()

Since 31 January 2025, total accumulated rainfall has broken multiple records across north Queensland. Record-crushing level of rain was observed in the coastal area north from Townsville.