Idealized Case 1: Tracking of a Test Blob in 2D

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

  1. Input Data

  2. Feature Detection

  3. Tracking / Trajectory Linking

  4. Segmentation

  5. 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.5.3

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

[2]:

import matplotlib.pyplot as plt

%matplotlib inline

import numpy as np

For a better readability of the graphs:

[3]:

import seaborn as sns

sns.set_context("talk")

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

1. Input Data

There are several utilities implemented in tobac to create simple examples of such arrays. In this tutorial we will use the function make_simple_sample_data_2D() to create a moving test blob in 2D:

[4]:

test_data = tobac.testing.make_simple_sample_data_2D(data_type="xarray")
test_data

[4]:

<xarray.DataArray 'w' (time: 100, y: 50, x: 100)> Size: 4MB
[500000 values with dtype=float64]
Coordinates:
  * time       (time) datetime64[ns] 800B 2000-01-01T12:00:00 ... 2000-01-01T...
  * y          (y) float64 400B 0.0 1e+03 2e+03 ... 4.7e+04 4.8e+04 4.9e+04
  * x          (x) float64 800B 0.0 1e+03 2e+03 ... 9.7e+04 9.8e+04 9.9e+04
    latitude   (y, x) float64 40kB ...
    longitude  (y, x) float64 40kB ...
Attributes:
    units:    m s-1

As you can see our generated data describes a field called ‘w’ with the unit m/s at 100, 50 and 100 datapoints of time, x and y. Additionally, the data contains the latitude and longitude coordinates of the field values. To access the values of ‘w’ in the first timeframe, we can use

[5]:

test_data.data[0]

[5]:

array([[3.67879441e+00, 4.04541885e+00, 4.40431655e+00, ...,
        2.22319774e-16, 9.26766698e-17, 3.82489752e-17],
       [4.04541885e+00, 4.44858066e+00, 4.84324569e+00, ...,
        2.44475908e-16, 1.01912721e-16, 4.20608242e-17],
       [4.40431655e+00, 4.84324569e+00, 5.27292424e+00, ...,
        2.66165093e-16, 1.10954118e-16, 4.57923372e-17],
       ...,
       [6.45813368e-03, 7.10174389e-03, 7.73178977e-03, ...,
        3.90282972e-19, 1.62694148e-19, 6.71461808e-20],
       [4.43860604e-03, 4.88095244e-03, 5.31397622e-03, ...,
        2.68237303e-19, 1.11817944e-19, 4.61488502e-20],
       [3.02025230e-03, 3.32124719e-03, 3.61589850e-03, ...,
        1.82522243e-19, 7.60865909e-20, 3.14020145e-20]])

which is then just an array of numbers of the described shape.

To visualize the data, we can plot individual time frames using matplotlib per imshow:

[6]:

fig, axs = plt.subplots(ncols=1, nrows=3, figsize=(12, 16), sharey=True)
plt.subplots_adjust(hspace=0.5)

for i, itime in enumerate([0, 20, 40]):
    # plot the 2D blob field in colors
    test_data.isel(time=itime).plot(ax=axs[i])

    axs[i].set_title(f"timeframe = {itime}")
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_15_0.png

This tells as that our data is a single moving blob, which is what we are going the analyze with tobac now.

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:

[7]:

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:

[8]:

test_data.max()

[8]:

<xarray.DataArray 'w' ()> Size: 8B
array(10.)

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.

[9]:

threshold = 9

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

[10]:

%%capture
features = tobac.feature_detection_multithreshold(test_data, dxy, threshold)

Let’s inspect the resulting object:

[11]:

features

[11]:
frame idx hdim_1 hdim_2 num threshold_value feature time timestr projection_y_coordinate projection_x_coordinate latitude longitude
0 0 1 10.000000 10.000000 69 9 1 2000-01-01 12:00:00 2000-01-01 12:00:00 10000.000000 10000.000000 24.100000 150.100000
1 1 1 10.939394 11.848485 66 9 2 2000-01-01 12:01:00 2000-01-01 12:01:00 10939.393939 11848.484848 24.118485 150.109394
2 2 1 11.707692 13.661538 65 9 3 2000-01-01 12:02:00 2000-01-01 12:02:00 11707.692308 13661.538462 24.136615 150.117077
3 3 1 12.569231 15.353846 65 9 4 2000-01-01 12:03:00 2000-01-01 12:03:00 12569.230769 15353.846154 24.153538 150.125692
4 4 1 13.200000 17.107692 65 9 5 2000-01-01 12:04:00 2000-01-01 12:04:00 13200.000000 17107.692308 24.171077 150.132000
5 5 1 14.184615 19.000000 65 9 6 2000-01-01 12:05:00 2000-01-01 12:05:00 14184.615385 19000.000000 24.190000 150.141846
6 6 1 15.060606 20.848485 66 9 7 2000-01-01 12:06:00 2000-01-01 12:06:00 15060.606061 20848.484848 24.208485 150.150606
7 7 1 15.953125 22.765625 64 9 8 2000-01-01 12:07:00 2000-01-01 12:07:00 15953.125000 22765.625000 24.227656 150.159531
8 8 1 16.707692 24.338462 65 9 9 2000-01-01 12:08:00 2000-01-01 12:08:00 16707.692308 24338.461538 24.243385 150.167077
9 9 1 17.676923 26.169231 65 9 10 2000-01-01 12:09:00 2000-01-01 12:09:00 17676.923077 26169.230769 24.261692 150.176769
10 10 1 18.184615 28.000000 65 9 11 2000-01-01 12:10:00 2000-01-01 12:10:00 18184.615385 28000.000000 24.280000 150.181846
11 11 1 19.171875 29.828125 64 9 12 2000-01-01 12:11:00 2000-01-01 12:11:00 19171.875000 29828.125000 24.298281 150.191719
12 12 1 20.046875 31.765625 64 9 13 2000-01-01 12:12:00 2000-01-01 12:12:00 20046.875000 31765.625000 24.317656 150.200469
13 13 1 20.953125 33.234375 64 9 14 2000-01-01 12:13:00 2000-01-01 12:13:00 20953.125000 33234.375000 24.332344 150.209531
14 14 1 21.828125 35.171875 64 9 15 2000-01-01 12:14:00 2000-01-01 12:14:00 21828.125000 35171.875000 24.351719 150.218281
15 15 1 22.815385 37.000000 65 9 16 2000-01-01 12:15:00 2000-01-01 12:15:00 22815.384615 37000.000000 24.370000 150.228154
16 16 1 23.323077 38.830769 65 9 17 2000-01-01 12:16:00 2000-01-01 12:16:00 23323.076923 38830.769231 24.388308 150.233231
17 17 1 24.292308 40.661538 65 9 18 2000-01-01 12:17:00 2000-01-01 12:17:00 24292.307692 40661.538462 24.406615 150.242923
18 18 1 25.046875 42.234375 64 9 19 2000-01-01 12:18:00 2000-01-01 12:18:00 25046.875000 42234.375000 24.422344 150.250469
19 19 1 25.939394 44.151515 66 9 20 2000-01-01 12:19:00 2000-01-01 12:19:00 25939.393939 44151.515152 24.441515 150.259394
20 20 1 26.815385 46.000000 65 9 21 2000-01-01 12:20:00 2000-01-01 12:20:00 26815.384615 46000.000000 24.460000 150.268154
21 21 1 27.800000 47.892308 65 9 22 2000-01-01 12:21:00 2000-01-01 12:21:00 27800.000000 47892.307692 24.478923 150.278000
22 22 1 28.430769 49.646154 65 9 23 2000-01-01 12:22:00 2000-01-01 12:22:00 28430.769231 49646.153846 24.496462 150.284308
23 23 1 29.292308 51.338462 65 9 24 2000-01-01 12:23:00 2000-01-01 12:23:00 29292.307692 51338.461538 24.513385 150.292923
24 24 1 30.060606 53.151515 66 9 25 2000-01-01 12:24:00 2000-01-01 12:24:00 30060.606061 53151.515152 24.531515 150.300606
25 25 1 31.000000 55.000000 69 9 26 2000-01-01 12:25:00 2000-01-01 12:25:00 31000.000000 55000.000000 24.550000 150.310000
26 26 1 31.939394 56.848485 66 9 27 2000-01-01 12:26:00 2000-01-01 12:26:00 31939.393939 56848.484848 24.568485 150.319394
27 27 1 32.707692 58.661538 65 9 28 2000-01-01 12:27:00 2000-01-01 12:27:00 32707.692308 58661.538462 24.586615 150.327077
28 28 1 33.569231 60.353846 65 9 29 2000-01-01 12:28:00 2000-01-01 12:28:00 33569.230769 60353.846154 24.603538 150.335692
29 29 1 34.200000 62.107692 65 9 30 2000-01-01 12:29:00 2000-01-01 12:29:00 34200.000000 62107.692308 24.621077 150.342000
30 30 1 35.184615 64.000000 65 9 31 2000-01-01 12:30:00 2000-01-01 12:30:00 35184.615385 64000.000000 24.640000 150.351846
31 31 1 36.060606 65.848485 66 9 32 2000-01-01 12:31:00 2000-01-01 12:31:00 36060.606061 65848.484848 24.658485 150.360606
32 32 1 36.953125 67.765625 64 9 33 2000-01-01 12:32:00 2000-01-01 12:32:00 36953.125000 67765.625000 24.677656 150.369531
33 33 1 37.707692 69.338462 65 9 34 2000-01-01 12:33:00 2000-01-01 12:33:00 37707.692308 69338.461538 24.693385 150.377077
34 34 1 38.676923 71.169231 65 9 35 2000-01-01 12:34:00 2000-01-01 12:34:00 38676.923077 71169.230769 24.711692 150.386769
35 35 1 39.184615 73.000000 65 9 36 2000-01-01 12:35:00 2000-01-01 12:35:00 39184.615385 73000.000000 24.730000 150.391846
36 36 1 40.171875 74.828125 64 9 37 2000-01-01 12:36:00 2000-01-01 12:36:00 40171.875000 74828.125000 24.748281 150.401719
37 37 1 41.046875 76.765625 64 9 38 2000-01-01 12:37:00 2000-01-01 12:37:00 41046.875000 76765.625000 24.767656 150.410469
38 38 1 41.953125 78.234375 64 9 39 2000-01-01 12:38:00 2000-01-01 12:38:00 41953.125000 78234.375000 24.782344 150.419531
39 39 1 42.828125 80.171875 64 9 40 2000-01-01 12:39:00 2000-01-01 12:39:00 42828.125000 80171.875000 24.801719 150.428281
40 40 1 43.815385 82.000000 65 9 41 2000-01-01 12:40:00 2000-01-01 12:40:00 43815.384615 82000.000000 24.820000 150.438154
41 41 1 44.462687 83.820896 67 9 42 2000-01-01 12:41:00 2000-01-01 12:41:00 44462.686567 83820.895522 24.838209 150.444627
42 42 1 45.292308 85.661538 65 9 43 2000-01-01 12:42:00 2000-01-01 12:42:00 45292.307692 85661.538462 24.856615 150.452923
43 43 1 45.836066 87.278689 61 9 44 2000-01-01 12:43:00 2000-01-01 12:43:00 45836.065574 87278.688525 24.872787 150.458361
44 44 1 46.254545 89.145455 55 9 45 2000-01-01 12:44:00 2000-01-01 12:44:00 46254.545455 89145.454545 24.891455 150.462545
45 45 1 46.770833 91.000000 48 9 46 2000-01-01 12:45:00 2000-01-01 12:45:00 46770.833333 91000.000000 24.910000 150.467708
46 46 1 47.256410 93.000000 39 9 47 2000-01-01 12:46:00 2000-01-01 12:46:00 47256.410256 93000.000000 24.930000 150.472564
47 47 1 47.500000 94.764706 34 9 48 2000-01-01 12:47:00 2000-01-01 12:47:00 47500.000000 94764.705882 24.947647 150.475000
48 48 1 47.782609 96.130435 23 9 49 2000-01-01 12:48:00 2000-01-01 12:48:00 47782.608696 96130.434783 24.961304 150.477826
49 49 1 48.153846 97.230769 13 9 50 2000-01-01 12:49:00 2000-01-01 12:49:00 48153.846154 97230.769231 24.972308 150.481538
50 50 1 48.600000 98.200000 5 9 51 2000-01-01 12:50:00 2000-01-01 12:50:00 48600.000000 98200.000000 24.982000 150.486000

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

We can plot the detected feature positions on top of the colored blob.

[12]:

fig, axs = plt.subplots(ncols=1, nrows=3, figsize=(12, 16), sharey=True)
plt.subplots_adjust(hspace=0.5)

for i, itime in enumerate([0, 20, 40]):
    # plot the 2D blob field in colors
    test_data.isel(time=itime).plot(ax=axs[i])

    # plot the detected feature as black cross
    f = features.loc[[itime]]
    f.plot.scatter(
        x="projection_x_coordinate",
        y="projection_y_coordinate",
        s=40,
        ax=axs[i],
        color="black",
        marker="x",
    )

    axs[i].set_title(f"timeframe = {itime}")
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_30_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.

[13]:

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

Frame 50: 1 trajectories present.

fails without v_max

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

[14]:

trajectories

[14]:
frame idx hdim_1 hdim_2 num threshold_value feature time timestr projection_y_coordinate projection_x_coordinate latitude longitude cell time_cell
0 0 1 10.000000 10.000000 69 9 1 2000-01-01 12:00:00 2000-01-01 12:00:00 10000.000000 10000.000000 24.100000 150.100000 1 0 days 00:00:00
1 1 1 10.939394 11.848485 66 9 2 2000-01-01 12:01:00 2000-01-01 12:01:00 10939.393939 11848.484848 24.118485 150.109394 1 0 days 00:01:00
2 2 1 11.707692 13.661538 65 9 3 2000-01-01 12:02:00 2000-01-01 12:02:00 11707.692308 13661.538462 24.136615 150.117077 1 0 days 00:02:00
3 3 1 12.569231 15.353846 65 9 4 2000-01-01 12:03:00 2000-01-01 12:03:00 12569.230769 15353.846154 24.153538 150.125692 1 0 days 00:03:00
4 4 1 13.200000 17.107692 65 9 5 2000-01-01 12:04:00 2000-01-01 12:04:00 13200.000000 17107.692308 24.171077 150.132000 1 0 days 00:04:00
5 5 1 14.184615 19.000000 65 9 6 2000-01-01 12:05:00 2000-01-01 12:05:00 14184.615385 19000.000000 24.190000 150.141846 1 0 days 00:05:00
6 6 1 15.060606 20.848485 66 9 7 2000-01-01 12:06:00 2000-01-01 12:06:00 15060.606061 20848.484848 24.208485 150.150606 1 0 days 00:06:00
7 7 1 15.953125 22.765625 64 9 8 2000-01-01 12:07:00 2000-01-01 12:07:00 15953.125000 22765.625000 24.227656 150.159531 1 0 days 00:07:00
8 8 1 16.707692 24.338462 65 9 9 2000-01-01 12:08:00 2000-01-01 12:08:00 16707.692308 24338.461538 24.243385 150.167077 1 0 days 00:08:00
9 9 1 17.676923 26.169231 65 9 10 2000-01-01 12:09:00 2000-01-01 12:09:00 17676.923077 26169.230769 24.261692 150.176769 1 0 days 00:09:00
10 10 1 18.184615 28.000000 65 9 11 2000-01-01 12:10:00 2000-01-01 12:10:00 18184.615385 28000.000000 24.280000 150.181846 1 0 days 00:10:00
11 11 1 19.171875 29.828125 64 9 12 2000-01-01 12:11:00 2000-01-01 12:11:00 19171.875000 29828.125000 24.298281 150.191719 1 0 days 00:11:00
12 12 1 20.046875 31.765625 64 9 13 2000-01-01 12:12:00 2000-01-01 12:12:00 20046.875000 31765.625000 24.317656 150.200469 1 0 days 00:12:00
13 13 1 20.953125 33.234375 64 9 14 2000-01-01 12:13:00 2000-01-01 12:13:00 20953.125000 33234.375000 24.332344 150.209531 1 0 days 00:13:00
14 14 1 21.828125 35.171875 64 9 15 2000-01-01 12:14:00 2000-01-01 12:14:00 21828.125000 35171.875000 24.351719 150.218281 1 0 days 00:14:00
15 15 1 22.815385 37.000000 65 9 16 2000-01-01 12:15:00 2000-01-01 12:15:00 22815.384615 37000.000000 24.370000 150.228154 1 0 days 00:15:00
16 16 1 23.323077 38.830769 65 9 17 2000-01-01 12:16:00 2000-01-01 12:16:00 23323.076923 38830.769231 24.388308 150.233231 1 0 days 00:16:00
17 17 1 24.292308 40.661538 65 9 18 2000-01-01 12:17:00 2000-01-01 12:17:00 24292.307692 40661.538462 24.406615 150.242923 1 0 days 00:17:00
18 18 1 25.046875 42.234375 64 9 19 2000-01-01 12:18:00 2000-01-01 12:18:00 25046.875000 42234.375000 24.422344 150.250469 1 0 days 00:18:00
19 19 1 25.939394 44.151515 66 9 20 2000-01-01 12:19:00 2000-01-01 12:19:00 25939.393939 44151.515152 24.441515 150.259394 1 0 days 00:19:00
20 20 1 26.815385 46.000000 65 9 21 2000-01-01 12:20:00 2000-01-01 12:20:00 26815.384615 46000.000000 24.460000 150.268154 1 0 days 00:20:00
21 21 1 27.800000 47.892308 65 9 22 2000-01-01 12:21:00 2000-01-01 12:21:00 27800.000000 47892.307692 24.478923 150.278000 1 0 days 00:21:00
22 22 1 28.430769 49.646154 65 9 23 2000-01-01 12:22:00 2000-01-01 12:22:00 28430.769231 49646.153846 24.496462 150.284308 1 0 days 00:22:00
23 23 1 29.292308 51.338462 65 9 24 2000-01-01 12:23:00 2000-01-01 12:23:00 29292.307692 51338.461538 24.513385 150.292923 1 0 days 00:23:00
24 24 1 30.060606 53.151515 66 9 25 2000-01-01 12:24:00 2000-01-01 12:24:00 30060.606061 53151.515152 24.531515 150.300606 1 0 days 00:24:00
25 25 1 31.000000 55.000000 69 9 26 2000-01-01 12:25:00 2000-01-01 12:25:00 31000.000000 55000.000000 24.550000 150.310000 1 0 days 00:25:00
26 26 1 31.939394 56.848485 66 9 27 2000-01-01 12:26:00 2000-01-01 12:26:00 31939.393939 56848.484848 24.568485 150.319394 1 0 days 00:26:00
27 27 1 32.707692 58.661538 65 9 28 2000-01-01 12:27:00 2000-01-01 12:27:00 32707.692308 58661.538462 24.586615 150.327077 1 0 days 00:27:00
28 28 1 33.569231 60.353846 65 9 29 2000-01-01 12:28:00 2000-01-01 12:28:00 33569.230769 60353.846154 24.603538 150.335692 1 0 days 00:28:00
29 29 1 34.200000 62.107692 65 9 30 2000-01-01 12:29:00 2000-01-01 12:29:00 34200.000000 62107.692308 24.621077 150.342000 1 0 days 00:29:00
30 30 1 35.184615 64.000000 65 9 31 2000-01-01 12:30:00 2000-01-01 12:30:00 35184.615385 64000.000000 24.640000 150.351846 1 0 days 00:30:00
31 31 1 36.060606 65.848485 66 9 32 2000-01-01 12:31:00 2000-01-01 12:31:00 36060.606061 65848.484848 24.658485 150.360606 1 0 days 00:31:00
32 32 1 36.953125 67.765625 64 9 33 2000-01-01 12:32:00 2000-01-01 12:32:00 36953.125000 67765.625000 24.677656 150.369531 1 0 days 00:32:00
33 33 1 37.707692 69.338462 65 9 34 2000-01-01 12:33:00 2000-01-01 12:33:00 37707.692308 69338.461538 24.693385 150.377077 1 0 days 00:33:00
34 34 1 38.676923 71.169231 65 9 35 2000-01-01 12:34:00 2000-01-01 12:34:00 38676.923077 71169.230769 24.711692 150.386769 1 0 days 00:34:00
35 35 1 39.184615 73.000000 65 9 36 2000-01-01 12:35:00 2000-01-01 12:35:00 39184.615385 73000.000000 24.730000 150.391846 1 0 days 00:35:00
36 36 1 40.171875 74.828125 64 9 37 2000-01-01 12:36:00 2000-01-01 12:36:00 40171.875000 74828.125000 24.748281 150.401719 1 0 days 00:36:00
37 37 1 41.046875 76.765625 64 9 38 2000-01-01 12:37:00 2000-01-01 12:37:00 41046.875000 76765.625000 24.767656 150.410469 1 0 days 00:37:00
38 38 1 41.953125 78.234375 64 9 39 2000-01-01 12:38:00 2000-01-01 12:38:00 41953.125000 78234.375000 24.782344 150.419531 1 0 days 00:38:00
39 39 1 42.828125 80.171875 64 9 40 2000-01-01 12:39:00 2000-01-01 12:39:00 42828.125000 80171.875000 24.801719 150.428281 1 0 days 00:39:00
40 40 1 43.815385 82.000000 65 9 41 2000-01-01 12:40:00 2000-01-01 12:40:00 43815.384615 82000.000000 24.820000 150.438154 1 0 days 00:40:00
41 41 1 44.462687 83.820896 67 9 42 2000-01-01 12:41:00 2000-01-01 12:41:00 44462.686567 83820.895522 24.838209 150.444627 1 0 days 00:41:00
42 42 1 45.292308 85.661538 65 9 43 2000-01-01 12:42:00 2000-01-01 12:42:00 45292.307692 85661.538462 24.856615 150.452923 1 0 days 00:42:00
43 43 1 45.836066 87.278689 61 9 44 2000-01-01 12:43:00 2000-01-01 12:43:00 45836.065574 87278.688525 24.872787 150.458361 1 0 days 00:43:00
44 44 1 46.254545 89.145455 55 9 45 2000-01-01 12:44:00 2000-01-01 12:44:00 46254.545455 89145.454545 24.891455 150.462545 1 0 days 00:44:00
45 45 1 46.770833 91.000000 48 9 46 2000-01-01 12:45:00 2000-01-01 12:45:00 46770.833333 91000.000000 24.910000 150.467708 1 0 days 00:45:00
46 46 1 47.256410 93.000000 39 9 47 2000-01-01 12:46:00 2000-01-01 12:46:00 47256.410256 93000.000000 24.930000 150.472564 1 0 days 00:46:00
47 47 1 47.500000 94.764706 34 9 48 2000-01-01 12:47:00 2000-01-01 12:47:00 47500.000000 94764.705882 24.947647 150.475000 1 0 days 00:47:00
48 48 1 47.782609 96.130435 23 9 49 2000-01-01 12:48:00 2000-01-01 12:48:00 47782.608696 96130.434783 24.961304 150.477826 1 0 days 00:48:00
49 49 1 48.153846 97.230769 13 9 50 2000-01-01 12:49:00 2000-01-01 12:49:00 48153.846154 97230.769231 24.972308 150.481538 1 0 days 00:49:00
50 50 1 48.600000 98.200000 5 9 51 2000-01-01 12:50:00 2000-01-01 12:50:00 48600.000000 98200.000000 24.982000 150.486000 1 0 days 00:50: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:

[15]:

track_mask = trajectories["cell"] == 1.0

This mask can then be used - to select the track:

[16]:

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

  • and to show the track in our plots:

[17]:

fig, axs = plt.subplots(ncols=1, nrows=3, figsize=(12, 16), sharey=True)
plt.subplots_adjust(hspace=0.5)
# fig.tight_layout()


for i, itime in enumerate([0, 20, 40]):
    # plot the 2D blob field in colors
    test_data.isel(time=itime).plot(ax=axs[i])

    # plot the track
    track.plot.scatter(
        x="projection_x_coordinate",
        y="projection_y_coordinate",
        ax=axs[i],
        color="red",
        marker="o",
        alpha=0.2,
    )

    # plot the detected feature as black cross
    f = features.loc[[itime]]
    f.plot.scatter(
        x="projection_x_coordinate",
        y="projection_y_coordinate",
        s=40,
        ax=axs[i],
        color="black",
        marker="x",
    )

    axs[i].set_title(f"timeframe = {itime}")
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_42_0.png

4. Segmentation

Another step after the feature detection and tracking is the segmentation of the data. This means we find the area surrounding our features belonging to the same cluster. Logically, we now need the already detected features as an additional input.

Just a small thought here:

We often introduce the different steps as 1.feature detection, 2.segmentation and 3. tracking, so having the segmenation step actually before the trajectory linking. The order actually does not matter since the current segmentation approach is based on the feature detection output only and can be done before or after the linking.

On further extensions:

  • Also, one could actually use the output of the segmentation (segments Dataset in the example below) as input for the tracking with the advantage that information on the area (ncells) is added in the dataframe.

  • One could also use the output of the tracking in the segmentation (trajectories Dataset from above) with the advantage that mask will contain only the features that are also linked to trajectories.

These possibilities exist and could be utlized by you …

[18]:

%%capture
segment_labels, segments = tobac.segmentation_2D(features, test_data, dxy, threshold=9)

As the name implies, the first object returned is an array in which the segmented areas belonging to each feature have the same label value as that feature. The second output is a dataframe of the detected features updated with information about their segmented regions (currently only the number of pixels segmented)

[19]:

segment_labels

[19]:

<xarray.DataArray 'segmentation_mask' (time: 100, y: 50, x: 100)> Size: 2MB
[500000 values with dtype=int32]
Coordinates:
  * time       (time) datetime64[ns] 800B 2000-01-01T12:00:00 ... 2000-01-01T...
  * y          (y) float64 400B 0.0 1e+03 2e+03 ... 4.7e+04 4.8e+04 4.9e+04
  * x          (x) float64 800B 0.0 1e+03 2e+03 ... 9.7e+04 9.8e+04 9.9e+04
    latitude   (y, x) float64 40kB ...
    longitude  (y, x) float64 40kB ...
Attributes:
    long_name:  segmentation_mask

Since True/False corresponds to 1/0 in python, we can visualize the segments with a simple contour plot:

[20]:

fig, axs = plt.subplots(ncols=1, nrows=3, figsize=(12, 16), sharey=True)
plt.subplots_adjust(hspace=0.5)


for i, itime in enumerate([0, 20, 40]):
    # plot the 2D blob field in colors
    test_data.isel(time=itime).plot(ax=axs[i])

    # plot the mask outline
    segment_labels.isel(time=itime).plot.contour(levels=[0.5], ax=axs[i], colors="k")

    # plot the detected feature as black cross
    f = features.loc[[itime]]
    f.plot.scatter(
        x="projection_x_coordinate",
        y="projection_y_coordinate",
        s=40,
        ax=axs[i],
        color="black",
        marker="x",
    )

    axs[i].set_title(f"timeframe = {itime}")
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_48_0.png

Keep in mind that the area of the resulting segments crucially depends on the defined thresholds.

5. Statistical Analysis

Blob Velocity / Motion Speed

Now different functions of tobacs analysis module can be used to calculate or plot properties of the tracks and the features. For example the velocity of the feature along the track can be calulated via:

[21]:

vel = tobac.calculate_velocity(track)

The expected velocity (hard-coded in the test function) is:

[22]:

expected_velocity = np.sqrt(30**2 + 14**2)

Plotting the velocity vs the timeframe will give us

[23]:

plt.figure(figsize=(10, 5))
plt.tight_layout()
plt.plot(vel["frame"], vel["v"])

plt.xlabel("timeframe")
plt.ylabel("velocity [m/s]")
plt.grid()

plt.axhline(expected_velocity, color="darkgreen", lw=5, alpha=0.5)
sns.despine()
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_57_0.png

We can also create an histogram of the detected velocities:

[24]:

hist, edges = tobac.velocity_histogram(
    track,
    bin_edges=np.arange(14, 43, 3),
)

[25]:

width = 0.9 * (edges[1] - edges[0])
center = (edges[:-1] + edges[1:]) / 2

plt.tight_layout()
plt.bar(center, hist, width=width)
plt.ylabel("count")
plt.xlabel("velocity [m/s]")

plt.axvline(expected_velocity, color="darkgreen", lw=5, alpha=0.5)
sns.despine()
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_60_0.png

Area

The area of the features can also be calculated and plotted throughout time:

[26]:

area = tobac.calculate_area(features, segment_labels)

[27]:

blob_magitude = 10  # also hard-code in the test
blob_sigma = 10e3

normalized_circle_radius = np.sqrt(np.log(blob_magitude / threshold))
absolute_circle_radius = np.sqrt(2) * blob_sigma * normalized_circle_radius
expected_area = np.pi * absolute_circle_radius**2

[28]:

plt.figure(figsize=(10, 5))
plt.tight_layout()
plt.plot(area["frame"], area["area"])
plt.xlabel("timeframe")
plt.ylabel(r"area [$m^2$]")
plt.grid()

plt.axhline(expected_area, color="darkgreen", lw=5, alpha=0.5)
sns.despine()
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_64_0.png

Lifetime

Another interesting property for real data are the lifetimes of our features. Tobac can also produce a histogram of this:

[29]:

hist, bins, centers = tobac.lifetime_histogram(
    track, bin_edges=np.arange(0, 200, 20)
)

[30]:

plt.tight_layout()
plt.bar(centers, hist, width=width)
plt.ylabel("count")
plt.xlabel("lifetime")
plt.show()
../../_images/examples4doc_Basics_Idealized-Case-1_Tracking-of-a-Test-Blob-in-2D_67_0.png

We can deduce, that our singular feature in the data has a lifetime of 50.