Batching Tutorial

This notebook works through a TCTrack workflow using the batching utility to run the Tempest Extremes tracker over the ERA5 reanalysis example data in batches, rather than all in one go.

Batching splits the time range into a series of batches, running the tracker over each in turn. This is useful for long datasets, where processing the data all at once may use too much memory or disk space, or when the data for each period needs to be retrieved separately.

We start by installing the tracker, then obtain and preprocess the example data, run the tracking in batches, and finally visualise the resulting tracks.

Setup

This notebook assumes that you have already run the Tempest Extremes notebook to install the tracking algorithm software and retrieve the data. In which case, we only need to update the PATH.

import os

os.environ["PATH"] = f"{os.getcwd()}/tempestextremes/build/bin:{os.environ['PATH']}"

Batching Function

The full documentation for the batching function can be found here.

For each batch period it will select the data for the appropriate time period and save it to a batch_[i]/ directory. Any preprocessing is then performed before running the tracker (using run_tracker). After the trajectories have been found for all batches it will combine them into a single file.

Briefly, the required input arguments are:

• A tracker object.

• A list of input files.

• An interval for each batch ("month", "year", or a cf.TimeDuration object).

• The start and end dates in "YYYY-MM-DD".

And the additional, optional arguments are:

• A list of preprocessing functions to perform on the input data.

• A function to manually retrieve data for each batch.

• A list of input files to pass to the tracker (after preprocessing).

• Additional configuration options.

For this tutorial we will use the seven-day data with a batch interval of four days. i.e. two batch periods, the first is the full four days, the second is the remaining three days. Typically, much longer batch intervals would be used.

First, we need to create the tracker. We will use the default parameter set.

import tctrack.tempest_extremes as te

detect_params, stitch_params = te.parameter_set("default", nc_names={"z": "z"})
tracker = te.TETracker(detect_params, stitch_params)

Input files

Now we can specify the list of input files.

Only the (non-surface) geopotential data is used directly without any preprocessing, so can be copied to the batch directory unchanged.

For the other files we can avoid unnecessarily saving them to the batch directory and instead keep them in memory by providing a dictionary with keys: "batch_file" (= None) and "store". The "store" values will be used to refer to the variables during the preprocessing.

Finally, the surface geopotential (used to find orography) is not time-varying and should not be subspaced to the batch time periods. To do this we use "time_varying" = False.

input_files = [
    "data/era5_z.nc",
    ("data/era5_sfc.nc", {"store": ["msl", "u10", "v10"], "batch_file": None}),
    ("data/era5_sfc_z.nc", {"store": "z", "batch_file": None, "time_varying": False}),
]

Preprocessing

For the preprocessing we will provide a list of functions to apply during each batch. Refer to the Tempest Extremes tutorial for details about why each function is applied.

For each preprocessing step, the second value is a dictionary of arguments to pass to the function. Note that "%BATCH%" can be used to refer to the batch directory.

The third value is another dictionary taking "use" to use the fields stored in memory, and "store" to keep the outputs in memory.

Here the preprocessing helper functions are used, but other functions can also be used. To be compatible with "use" / "store" these must accept as the first argument / return a cf.Field object (or a list of them).

import cf

from tctrack.preprocessing import (
    calculate_wind_speed,
    multiply_field,
    set_netcdf_info,
    squeeze_field,
)

GRAVITY = 9.80665  # Standard gravitational acceleration [m s-2]

preprocessing = [
    (cf.write, {"dataset_name": "%BATCH%/msl.nc"}, {"use": "msl"}),
    (
        calculate_wind_speed,
        {"nc_name": "si10", "output_file": "%BATCH%/si10.nc"},
        {"use": ["u10", "v10"]},
    ),
    (squeeze_field, {}, {"use": "z", "store": "z"}),
    (multiply_field, {"factor": 1 / GRAVITY}, {"use": "z", "store": "z"}),
    (
        set_netcdf_info,
        {
            "nc_name": "orog",
            "properties": {
                "standard_name": "surface_altitude",
                "long_name": "Surface Altitude",
                "units": "m",
            },
            "output_file": "%BATCH%/orog.nc",
        },
        {"use": "z"},
    ),
]

Batching config

We can use a dictionary to configure the behaviour of the batching function.

The default config values are used here, so it is not strictly necessary. However, it is included for clarity.

One option is to add a buffer period to prevent tracks from being cut off at the end of the batch period. Ideally this should be longer than most tropical cyclone tracks (so >= 10 days). If this is used, a short start buffer period prevents duplicate tracks (any tracks starting in this period are removed).

No buffer period is used here since the period is already so short.

from datetime import timedelta

config = {
    "output_dir": "tctrack_outputs",  # Where to save outputs
    "combine_outputs": True,  # Whether to produce a combined output file
    "delete_batch_dirs": False,  # Whether to delete the batch directories
    "buffer_period": None,  # The buffer period to include to avoid cut-off trajectories
    "start_buffer_period": timedelta(days=1),  # Not used if buffer_period is None
}

Run the batching function

Now we can run the batching function:

%%capture --no-stderr
from tctrack.utils import batching

batching(
    tracker,
    input_files,
    cf.TimeDuration(4, "days"),
    time_range=["1950-09-01", "1950-09-09"],
    preprocessing=preprocessing,
    config=config,
)

Arguments that are not used here are:

• retrieve_data: We already have the data for the full period, but this can be used to download the data for a batch at a time.

• tracker_inputs: By default it passes all of the files in the batch directory to the run_tracker function, which is what we wan.

Inspecting the output

The combined output file is a CF-NetCDF file in the output directory (tctrack_outputs/tracks.nc) containing the tracks from all the batches, along with variables such as wind speed and sea-level pressure along each track. The structure can be inspected with ncdump -h tctrack_outputs/tracks.nc, or read directly with cf-python:

import glob

print("Output files:", glob.glob("tctrack_outputs/tracks*.nc"))

tracks = cf.read("tctrack_outputs/tracks.nc")
print("Total tracks:", tracks[0].coordinate("trajectory").size)

Output files: ['tctrack_outputs/tracks_0.nc', 'tctrack_outputs/tracks.nc']
Total tracks: 1

As can be seen, only the first batch produced an output file and there was only one output track.

Compare this with the Tempest Extremes notebook where four tracks were identified. This is because the tracks are cut off at the edge of each batch period, and the resulting shorter tracks are filtered out.

We can test this by using a three-day buffer period to make the first batch cover the seven day period:

%%capture --no-stderr
batching(
    tracker,
    input_files,
    cf.TimeDuration(4, "days"),
    time_range=["1950-09-01", "1950-09-09"],
    preprocessing=preprocessing,
    config={"buffer_period": timedelta(days=3)},
)

print("Output files:", glob.glob("tctrack_outputs/tracks*.nc"))

tracks = cf.read("tctrack_outputs/tracks.nc")
print("Total tracks:", tracks[0].coordinate("trajectory").size)

Output files: ['tctrack_outputs/tracks_0.nc', 'tctrack_outputs/tracks.nc']
Total tracks: 4

Obviously, this makes the batching irrelevant in this case since it now runs the whole period in the first batch. However, it illustrates the concept for running the batching over longer time periods.

Visualising Results

Finally we can visualise the results by plotting the tracks on a map. These can be seen to match the results from running Tempest Extremes normally.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

# Extract the wind speed field and its coordinates as arrays
wind_field = tracks.select_field("wind_speed")
lats = wind_field.coordinate("latitude").array
lons = wind_field.coordinate("longitude").array
times = wind_field.coordinate("time").dtarray
intensity = wind_field.array

# Get intensity metadata for labels
intensity_name = wind_field.get_property("long_name")
intensity_units = wind_field.get_property("units")
min_intensity = np.nanmin(intensity)
max_intensity = np.nanmax(intensity)

plt.figure(figsize=(8, 4))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-180, 180, -90, 90])
ax.coastlines()
gl = ax.gridlines(draw_labels=True)
gl.top_labels = False
gl.right_labels = False
gl.left_labels = True
gl.bottom_labels = True

# Plot each trajectory
for i in range(lats.shape[0]):
    times_i = np.ma.compressed(times[i])
    label = (
        f"{times_i[0].strftime('%Y-%m-%d %H:%M')} to "
        f"{times_i[-1].strftime('%Y-%m-%d %H:%M')}"
    )
    mask = np.isfinite(lons[i])
    pl = ax.plot(
        lons[i, mask],
        lats[i, mask],
        "--",
        transform=ccrs.PlateCarree(),
        label=f"{label}",
    )
    sc = ax.scatter(
        lons[i],
        lats[i],
        c=intensity[i],
        cmap="viridis",
        s=40,
        vmin=min_intensity,
        vmax=max_intensity,
        transform=ccrs.PlateCarree(),
    )

plt.colorbar(sc, label=f"{intensity_name} ({intensity_units})")
plt.title(f"All Trajectories Colored by {intensity_name}")
plt.legend()
plt.tight_layout()
plt.show()