Scikit - Feature Tracking Cookbook

Overview¶

Within this notebook, we will cover:

How to filter continuous SLP data do identify features of interest.
How to use scikit-image’s label() function to identify coherent objects within the SLP field
How to track identified ETCs through successive time steps.

Prerequisites¶

Concepts	Importance	Notes
xarray	Necessary
matplotlib	Necessary
scikit-image	Useful	Mostly the the label function
Cartopy	Useful
IPython	Useful	For inline animation

Time to learn: 30 minutes

Import Packages¶

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
from skimage.measure import label, regionprops
import matplotlib.colors as mcolors
from matplotlib.cm import ScalarMappable
import cartopy.crs as ccrs

Let’s go through the basic steps to identify SLP features¶

We have a small subset of hourly ERA5 Reanalysis dataset for Sea Level Pressure for 2022 November 1.

ds = xr.open_dataset("data/SLP_ex.nc")
da = ds["SLP"]

Three values control what counts as a feature and how features are linked through time. We set them once here so the rest of the notebook stays readable.

threshold - ERA5 reports sea-level pressure in pascals, so this is 1005 hPa. Grid cells below this value are flagged as part of a low. Everything at or above it is background.
min_area- The smallest feature we keep, measured in grid cells. Below threshold blobs smaller than 20 cells are treated as noise and dropped. Since we haven’t smoothed the field yet, this is what stops single-cell specks from being counted as a feature.
min_overlap - The matching rule for tracking. A feature at one timestep is identified with a feature at the next when they share at least 20% overlap. Below that level, they’re treated as separate features. Lower values follow faster moving systems but risk false links and higher values require a feature to barely move between steps.

threshold = 100000
min_area = 20
min_overlap = 0.2

Get the latitude and longitude 2d meshgrid for plotting¶

To plot SLP with contour(), we first create 2D latitude–longitude grids. SLP : 1° × 1°

# ---- get lat/lon coords ----
# assumes dims like (time, lat, lon)
lat = da[da.dims[1]].values
lon = da[da.dims[2]].values

# if lat/lon are 1D, make 2D grids for contourf/contour
lon2d, lat2d = np.meshgrid(lon, lat)

Select low SLP features less than a threshold using scikit¶

The detect_features() function identifies low-pressure objects in a single SLP field. First, it creates a mask of all grid points where SLP is less than or equal to the chosen threshold. These grid points are treated as candidate low-pressure regions. The function then uses connected-component labeling to group neighboring grid points into separate features. Because connectivity=2 is used, grid cells that touch along edges or corners are considered connected. After the initial labeling, very small features are removed using the min_area limit, so noisy or insignificant patches are not included. Finally, the remaining features are labeled again to produce a clean label map, where 0 represents the background and each positive integer represents a separate low-pressure feature.

def detect_features(field, threshold=100000, min_area=20):
    mask = field <= threshold
    labels = label(mask, connectivity=2)

    for lab in np.unique(labels)[1:]:
        if np.sum(labels == lab) < min_area:
            labels[labels == lab] = 0

    labels = label(labels > 0, connectivity=2)
    return labels

Let’s look at the labelled feature plotted for one time step¶

Let’s take a look at the labelled feature at time = 0 which at 2022 November 1 00:00 UTC along with precipitation overlay.

t = 0
slp = da.isel(time=t).values

We call the detect_features() function with the parameters set above

labels = detect_features(slp, threshold=threshold, min_area=min_area)

levels = np.linspace(float(da.min()), float(da.max()), 16)

fig, ax = plt.subplots(
    figsize=(8, 4),
    subplot_kw={"projection": ccrs.PlateCarree()}
)

# colored SLP contour lines
ax.contour(
    lon2d, lat2d, slp,
    levels=levels,
    cmap="turbo",
    linewidths=0.9,
    transform=ccrs.PlateCarree()
)

ax.contour(
    lon2d,
    lat2d,
    labels.astype(int),
    levels=[0.5],
    colors="black",
    linewidths=1.6,
    transform=ccrs.PlateCarree()
)

for r in regionprops(labels):
    y, x = r.centroid      # row, column
    iy = int(round(y))
    ix = int(round(x))


    ax.plot(
        lon[ix],
        lat[iy],
        "ko",
        ms=3,
        transform=ccrs.PlateCarree()
    )

    ax.text(
        lon[ix],
        lat[iy],
        str(r.label),
        fontsize=8,
        color="black",
        ha="left",
        va="bottom",
        transform=ccrs.PlateCarree(),
        bbox=dict(facecolor="white", edgecolor="none", alpha=0.7, pad=1)
    )


ax.coastlines(color="0.55", linewidth=0.6)

gl = ax.gridlines(
    draw_labels=True,
    linewidth=0.4,
    color="0.5",
    alpha=0.25,
    linestyle="--"
)
gl.top_labels = False
gl.right_labels = False

#ax.set_title(f"Matplotlib contour-tracked SLP features | time index = {t}")

plt.show()

Let’s track SLP Features Over Time¶

We track the labelled SLP features from one time step to the next. Each detected feature is converted into a mask and compared with the masks from the previous time step. If a current feature overlaps strongly enough with a previous feature, it is treated as the same object and keeps the same track_id. If the overlap is too small, the feature is treated as a new object and is assigned a new track_id. For each time step, the code stores the feature ID, centroid location, and mask, which can later be used for plotting, animation, or further analysis of how the low-pressure systems move.

# ---- precompute tracked features ----
tracked = []
prev = {}
next_track_id = 1

for t in range(da.sizes["time"]):
    field = da.isel(time=t).values
    labels = detect_features(field, threshold=threshold, min_area=min_area)

    current = {}
    frame_info = []

    for r in regionprops(labels):
        mask = labels == r.label

        best_id = None
        best_overlap = 0.0

        for track_id, old_mask in prev.items():
            overlap = np.logical_and(mask, old_mask).sum() / mask.sum()
            if overlap > best_overlap:
                best_overlap = overlap
                best_id = track_id

        if best_overlap < min_overlap:
            best_id = next_track_id
            next_track_id += 1

        current[best_id] = mask
        frame_info.append({
            "track_id": best_id,
            "centroid_rc": r.centroid,   # row, col
            "mask": mask
        })

    tracked.append(frame_info)
    prev = current

Let’s look at the animation of the SLP features¶

# ---- color scale fixed across all frames ----
vmin = float(da.min())
vmax = float(da.max())
levels = np.linspace(vmin, vmax, 16)
norm = mcolors.Normalize(vmin=vmin, vmax=vmax)
cmap = "turbo"


# ---- animate ----
fig, ax = plt.subplots(figsize=(8, 4), subplot_kw={"projection": ccrs.PlateCarree()})

sm = ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
cbar = fig.colorbar(sm, ax=ax, pad=0.02)
cbar.set_label("SLP")


def update(t):
    ax.clear()

    field = da.isel(time=t).values

    # SLP as colored contour lines, not filled
    cs = ax.contour(
        lon2d, lat2d, field,
        levels=levels,
        cmap="turbo",          # or "rainbow", "viridis", "plasma"
        linewidths=0.9,
        transform=ccrs.PlateCarree()
    )

    # tracked feature outlines + labels
    for feat in tracked[t]:
        mask = feat["mask"]

        ax.contour(
            lon2d, lat2d, mask.astype(int),
            levels=[0.5],
            colors="black",
            linewidths=1.6,
            transform=ccrs.PlateCarree()
        )

        y, x = feat["centroid_rc"]
        iy = int(round(y))
        ix = int(round(x))

        ax.plot(
            lon[ix], lat[iy],
            "ko", ms=3,
            transform=ccrs.PlateCarree()
        )

        ax.text(
            lon[ix], lat[iy],
            str(feat["track_id"]),
            fontsize=8,
            color="black",
            ha="left",
            va="bottom",
            transform=ccrs.PlateCarree(),
            bbox=dict(facecolor="white", edgecolor="none", alpha=0.7, pad=1)
        )

    # continents/coastlines: faint, not black
    ax.coastlines(color="0.55", linewidth=0.6)

    # faint gridlines
    gl = ax.gridlines(
        draw_labels=True,
        linewidth=0.4,
        color="0.5",
        alpha=0.25,
        linestyle="--"
    )
    gl.top_labels = False
    gl.right_labels = False

    #ax.set_title(f"Scikit tracked SLP features | time index = {t}")

anim = FuncAnimation(
    fig,
    update,
    frames=da.sizes["time"],
    interval=300,
    blit=False,
)

plt.close(fig)
HTML(anim.to_jshtml())

This animation shows the evolution of SLP and precipitation through time using the overlap tracking method. The numeric annotations help follow individual systems across successive times, including cases where the detected object appears as nested within the larger low SLP features.

Summary¶

In this recipe, we identify low sea level pressure (SLP) features globally and track them through time. We use scikit-image.measure label() to detect coherent low-pressure objects. We track tracking identified features between time steps using spatial overlap percentage and made visualization for one time step and a hourly animation for duration of 1 day.

What’s next?¶

Adding smoothing or noise filtering before feature detection.
Comparing this contour-based method with scikit-image and scipy object-labeling approaches.
Extracting region properties such as area, centroid, perimeter, and intensity.
Testing sensitivity to different SLP thresholds.
Adding explicit split and merge detection during tracking.