Idealized Case 3: Tracking of a Test Blob in 3D#

This tutorial shows the most important steps of tracking with tobac using an idealized case:

Input Data
Feature Detection
Tracking / Trajectory Linking
Segmentation
Statistical Analysis

Import Libraries#

We start by importing tobac:

[1]:

import tobac

print("using tobac version", str(tobac.__version__))

# We add testing here to create test dataset (typically not needed in standard applications).
import tobac.testing

using tobac version 1.6.2

We will also need matplotlib in inline-mode for plotting and numpy:

[2]:

import matplotlib.pyplot as plt

%matplotlib inline

import numpy as np

Tobac works with a Python package called xarray, which introduces DataArrays. In a nutshell these are numpy-arrays with intuitive dimensions, coordinates and meta-data. For a more extensive description have a look at the xarray Documentation.

[3]:

import os
CI = os.environ.get("CI") == "true"

2. Feature Detection#

The first step of the general working routine of tobac is the identification of features. This essentially means finding the maxima or minima of the data.

To use the according functions of tobac we need to specify:

the thresholds below/above the features are detected
the spacing of our data

The spacing of the temporal and spatial dimension can be extracted from the data using a build-in utility:

[10]:

dxy, dt = tobac.get_spacings(test_data)

To get an idea of which order of magnitude our thresholds should be, we check the maximum of our data:

[11]:

test_data.max()

[11]:

<xarray.DataArray 'w' ()> Size: 8B
array(10.)
Attributes:
    units:    m s-1

Since we know that our data will only have one maximum it is reasoable to choose 9 as our threshold, but keep in mind that we could also add multiple values here if our data would be more complex.

[12]:

threshold = 9

Now we are ready to apply the feature detection algorithm. Notice that this is a minimal input. The function has several other options we will cover in later tutorials.

[13]:

# %%capture
features = tobac.feature_detection_multithreshold(test_data, dxy, threshold)
# This shows the first five rows
features.head(5)

[13]:

	frame	idx	vdim	hdim_1	hdim_2	num	threshold_value	feature	time	timestr	z	y	x	latitude	longitude
0	0	1	4.000000	10.000000	10.000000	33	9	1	2000-01-01 12:00:00	2000-01-01 12:00:00	4000.000000	10000.000000	10000.000000	24.100000	150.100000
1	1	1	5.840909	11.613636	13.545455	44	9	2	2000-01-01 12:02:00	2000-01-01 12:02:00	5840.909091	11613.636364	13545.454545	24.135455	150.116136
2	2	1	7.731707	13.390244	17.268293	41	9	3	2000-01-01 12:04:00	2000-01-01 12:04:00	7731.707317	13390.243902	17268.292683	24.172683	150.133902
3	3	1	9.658537	15.121951	20.756098	41	9	4	2000-01-01 12:06:00	2000-01-01 12:06:00	9658.536585	15121.951220	20756.097561	24.207561	150.151220
4	4	1	11.642857	16.761905	24.452381	42	9	5	2000-01-01 12:08:00	2000-01-01 12:08:00	11642.857143	16761.904762	24452.380952	24.244524	150.167619

The ouputs tells us that features were found in 25 frames (index 0 to 24) of our data. The variable idx is 1 for every frame, which means that only 1 feature was found in every time step, as we expected. vdim, hdim_1 and hdim_2 are the position of this feature with respect to the z, y and x-indices.

We can plot the detected feature positions in the 3D grid.

[14]:

# Create a PyVista Plotter object
pl = pv.Plotter()

# Create grid
x, y, z = np.meshgrid(test_data.x, test_data.y, test_data.z, indexing="ij")
grid = pv.StructuredGrid(x, y, z)

# Add the outline of the grid and axes
pl.add_mesh(grid.outline(), color="black", line_width=1)
pl.add_axes()
pl.show_bounds(**common_bounds_args)
pl.camera_position = camera_position

# Extract 3D coordinates of features
points = np.column_stack(
    [
        test_data.x.values[features["hdim_2"].astype(int)],
        test_data.y.values[features["hdim_1"].astype(int)],
        test_data.z.values[features["vdim"].astype(int)],
    ]
)

# Add grid outline and axes
pl.add_mesh(grid.outline(), color="black", line_width=1)
pl.add_axes()

labels = [f"({x}, {y}, {z})" for x, y, z in points]

# Collect visible points & labels for selected time steps
vis_points = points[time_steps]
vis_labels = [f"({x:.1f}, {y:.1f}, {z:.1f})" for x, y, z in vis_points]

# Add point labels for the selected time steps
pl.add_point_labels(vis_points, vis_labels, point_color="blue", point_size=10)

# Show the plot
if not CI:
    pl.show(auto_close=False)
pl.close()

../../_images/examples_Basics_Idealized-Case-3_Tracking-of-a-Test-Blob-in-3D_35_0.png

The function has succesfully detected the maximum of our data in the individual timeframes.

3. Trajectory Linking#

After we are done finding the features and associated segments for each frame it is necessary for further analysis to keep track of those elements troughout time. Linking is the tool for that. It connects the features of the timesteps which belong together. We are going to use the linking_trackpy() function here. The required inputs are the features, the two spacings and a maximum velocity of the features.

[15]:

trajectories = tobac.linking_trackpy(features, test_data, dt=dt, dxy=dxy, v_max=100)

Frame 24: 1 trajectories present.

Unsurprisingly, one trajectory was found. The returned object is another Dataset:

[16]:

# This shows the first five rows
trajectories.head(5)

[16]:

	frame	idx	vdim	hdim_1	hdim_2	num	threshold_value	feature	time	timestr	z	y	x	latitude	longitude	cell	time_cell
0	0	1	4.000000	10.000000	10.000000	33	9	1	2000-01-01 12:00:00	2000-01-01 12:00:00	4000.000000	10000.000000	10000.000000	24.100000	150.100000	1	0 days 00:00:00
1	1	1	5.840909	11.613636	13.545455	44	9	2	2000-01-01 12:02:00	2000-01-01 12:02:00	5840.909091	11613.636364	13545.454545	24.135455	150.116136	1	0 days 00:02:00
2	2	1	7.731707	13.390244	17.268293	41	9	3	2000-01-01 12:04:00	2000-01-01 12:04:00	7731.707317	13390.243902	17268.292683	24.172683	150.133902	1	0 days 00:04:00
3	3	1	9.658537	15.121951	20.756098	41	9	4	2000-01-01 12:06:00	2000-01-01 12:06:00	9658.536585	15121.951220	20756.097561	24.207561	150.151220	1	0 days 00:06:00
4	4	1	11.642857	16.761905	24.452381	42	9	5	2000-01-01 12:08:00	2000-01-01 12:08:00	11642.857143	16761.904762	24452.380952	24.244524	150.167619	1	0 days 00:08:00

The new variable cell now indexes the features in the different time steps. Therefore we can use it to create a mask for our moving feature:

[17]:

track_mask = trajectories["cell"] == 1.0

This mask can then be used

to select the track:

[18]:

track = trajectories.where(track_mask).dropna()

and to show the track in our 3D plots:

[19]:

# Create a PyVista Plotter object
pl = pv.Plotter()

# Create grid
x, y, z = np.meshgrid(test_data.x, test_data.y, test_data.z, indexing="ij")
grid = pv.StructuredGrid(x, y, z)

# Add grid outline and axes
pl.add_mesh(grid.outline(), color="black", line_width=1)
pl.add_axes()

# Extract 3D coordinates of features
points = np.column_stack(
    [
        test_data.x.values[features["hdim_2"].astype(int)],
        test_data.y.values[features["hdim_1"].astype(int)],
        test_data.z.values[features["vdim"].astype(int)],
    ]
)

tracks = np.column_stack(
    [
        test_data.x.values[track["hdim_2"].astype(int)],
        test_data.y.values[track["hdim_1"].astype(int)],
        test_data.z.values[track["vdim"].astype(int)],
    ]
)

labels = [f"({x}, {y}, {z})" for x, y, z in points]

# Collect visible points & labels for selected time steps
vis_points = points[time_steps]
vis_labels = [f"({x:.1f}, {y:.1f}, {z:.1f})" for x, y, z in vis_points]

# Add point labels for the selected time steps
pl.add_point_labels(vis_points, vis_labels, point_color="blue", point_size=10)

# Create the connectivity array for lines
lines = []
for i in range(len(tracks) - 1):
    lines.append([2, i, i + 1])

# Flatten the list of lines for the connectivity array
lines = np.array(lines).flatten()

# Create the PolyData object with the line segments
track_poly = pv.PolyData(tracks)
track_poly.lines = lines

# Add tracks to the plot
pl.add_mesh(track_poly, color="red", line_width=3)

volume_rendering(test_data, time_steps, pl)
pl.close()

../../_images/examples_Basics_Idealized-Case-3_Tracking-of-a-Test-Blob-in-3D_46_0.png

5. Bulk statistics for segmented features#

To characterise segmented objects, we can compute simple statistics of an underlying field over each labeled feature. This provides a direct link between the segmentation output and feature-level properties.

Using get_statistics_from_mask, we pass

the feature dataframe,
the segmentation mask, and
one or more fields (e.g. test_data),

and receive an updated dataframe with new columns containing the chosen bulk statistics (e.g. mean, total, max) for each feature.

[24]:

from tobac.utils.bulk_statistics import get_statistics_from_mask

# Start from the original features dataframe
features_input = features.copy()

# Segmentation mask and underlying 3D field from the previous steps
segmentation_mask = segment_labels
field_3d = test_data

# Define which statistics we want to compute for each feature
statistics = {}
statistics['mean'] = np.mean
statistics['total'] = np.sum
statistics['max'] = np.max

features_with_stats = get_statistics_from_mask(features_input, segmentation_mask, field_3d, statistic=statistics)

features_with_stats.head()

[24]:

	frame	idx	vdim	hdim_1	hdim_2	num	threshold_value	feature	time	timestr	z	y	x	latitude	longitude	mean	total	max
0	0	1	4.000000	10.000000	10.000000	33	9	1	2000-01-01 12:00:00	2000-01-01 12:00:00	4000.000000	10000.000000	10000.000000	24.100000	150.100000	9.333259	531.995777	10.0
1	1	1	5.840909	11.613636	13.545455	44	9	2	2000-01-01 12:02:00	2000-01-01 12:02:00	5840.909091	11613.636364	13545.454545	24.135455	150.116136	9.400239	470.011963	9.946384
2	2	1	7.731707	13.390244	17.268293	41	9	3	2000-01-01 12:04:00	2000-01-01 12:04:00	7731.707317	13390.243902	17268.292683	24.172683	150.133902	9.387958	478.785843	9.961036
3	3	1	9.658537	15.121951	20.756098	41	9	4	2000-01-01 12:06:00	2000-01-01 12:06:00	9658.536585	15121.951220	20756.097561	24.207561	150.151220	9.391814	469.590716	9.98018
4	4	1	11.642857	16.761905	24.452381	42	9	5	2000-01-01 12:08:00	2000-01-01 12:08:00	11642.857143	16761.904762	24452.380952	24.244524	150.167619	9.415543	451.946046	9.932072

The resulting dataframe features_with_stats contains the original feature information plus a set of new columns with the requested statistics (here: mean, total, and max) of test_data over the grid cells belonging to each segmented feature.

This bulk-statistics capability can be easily extended to multiple fields or different statistics, and it provides a flexible way to characterise the segmented objects.

[25]:

# For a simple overview, keep one row per feature ID
features_summary = (
    features_with_stats
    .sort_values("time")
    .drop_duplicates(subset=["feature"], keep="last")
)

plt.figure(figsize=(8, 5))
plt.scatter(
    features_summary["mean"],   # mean value of the field for each feature
    features_summary["max"],    # corresponding maximum
)

plt.xlabel("Mean value of test_data per feature")
plt.ylabel("Max value of test_data per feature")
plt.grid(True)
plt.tight_layout()

../../_images/examples_Basics_Idealized-Case-3_Tracking-of-a-Test-Blob-in-3D_58_0.png

Idealized Case 3: Tracking of a Test Blob in 3D#

Import Libraries#

1. Input Data#

2. Feature Detection#

3. Trajectory Linking#

4. Segmentation#

5. Bulk statistics for segmented features#

This Page