Gallery 6: Particle Tracking Plots¶
Plot particle tracking outputs from TUFLOW FV’s Particle Tracking Module.
This notebook is used in combination with the TUFLOW FV Python Toolbox (tfv) package.
These Python tools have been recently updated (May 2025) to use a familiar workflow to the hydrodynamic post processing tools with the .tfv accessor.
IMPORTANT: This module is under active development and has been released on the basis of obtaining early feedback. Please contact TUFLOW Support with issues / requests / feedback when using this new module.
from pathlib import Path
import matplotlib.pyplot as plt
import tfv.xarray
import xarray as xr
Load Results¶
We’re using two sets of results to show off different functionality.
folder = Path('../../data')
# Load up our first PTM and TFV Domain Result
ptm1 = xr.open_dataset(folder / 'PTM_002_ptm.nc', decode_timedelta=True).tfv
fv1 = xr.open_dataset(folder / 'PTM_002.nc').tfv
# Load up our second PTM and TFV Domain Result
ptm2 = xr.open_dataset(folder / 'PTM_005_ptm.nc', decode_timedelta=True).tfv
fv2 = xr.open_dataset(folder / 'PTM_005.nc').tfv
Scatter Plots¶
Plot Basic¶
The general purpose function that we use for filtering particles is .get_particles. That function is ran whenever you use a plotting function so you can just put in your filters here for ease of use.
# Pick a time to plot - we'll sue the same throughout the tutorial.
# You can enter a datetime e.g "2011-05-01 02:00";
# or just simple integer indexing, e.g., "-1" gives you the final timestep
time = -1
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Draw the mesh in the background
fv1.plot_mesh(ax=ax)
ptm1.plot(time=time, datum='height', limits=(0, 1), ax=ax, s=5, c='k', label='Bed (0m - 1m)')
ptm1.plot(time=time, datum='depth', limits=(0, 1), ax=ax, s=5, c='pink', label='Surface (0m - 1m)')
ax.legend()
plt.show()
Plot Groups¶
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Draw the mesh in the background
fv2.plot_mesh(ax=ax)
ptm2.plot(time=time, ax=ax, s=5, c='green', group_ids=1, label='Clay')
ptm2.plot(time=time, ax=ax, s=5, c='black', group_ids=2, label='Silt')
ptm2.plot(time=time, ax=ax, s=5, c='yellow', group_ids=3, label='Fine Sand')
ptm2.plot(time=time, ax=ax, s=5, c='brown', group_ids=4, label='Coarse Sand')
ax.legend()
plt.show()
Plot Groups State¶
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Draw the mesh in the background
fv1.plot_mesh(ax=ax)
ptm1.plot(time=time, ax=ax, s=5, c='blue', stats=1, label='In Suspension')
ptm1.plot(time=time, ax=ax, s=5, c='yellow', stats=2, label='Dry')
ptm1.plot(time=time, ax=ax, s=5, c='brown', stats=3, label='Bed')
ax.legend()
plt.show()
Plot Coloured¶
We can use the color_by argument to color the particles by a variable stored in the dataset. Examples include z, age, stats, group ids, etc.
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Draw the mesh in the background
fv1.plot_mesh(ax=ax)
ptm1.plot(time=time, ax=ax, s=5, color_by='z', cmap='jet')
plt.show()
Plot Highlight¶
Sometimes we just want to highlight a specific group or status. There’s a few ways we can do this, here’s a simple one if we know what group ids we have.
(You could totally just plot all particles and then overplot with the highlights too!)
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Draw the mesh in the background
fv1.plot_mesh(ax=ax)
ptm1.plot(time=time, ax=ax, s=5, c='k', group_ids=(1,3,4), label='Groups 1, 3 and 4')
ptm1.plot(time=time, ax=ax, s=10, c='y', group_ids=2, label='Group 2')
ax.legend()
plt.show()
Plot Highlight AND Status¶
Pretty much as above but let’s highlight a specific group and status. We’ll also turn off the mesh, but keep the boundary on.
fig, ax = plt.subplots()
ax.set_aspect('equal')
# Draw the mesh in the background
fv1.plot_mesh(ax=ax, mesh=False)
ptm1.plot(time=time, ax=ax, s=5, c='grey', alpha=0.25, label='All particles')
ptm1.plot(time=time, ax=ax, s=10, c='y', group_ids=2, stats=3, label='Special Particles')
ax.legend()
plt.show()
Grid Plots¶
The above plots are great when you have a managable set of particles, or you only want to show a single timestep.
However they quickly get bogged down when you have many many particles, such as order 100,000s to millions. In this case, grid plots are the way to go.
Pretty much all the above will still work, but you can supply plot_type=grid instead.
Plot Scatter vs Grid¶
fig, axes = plt.subplots(ncols=3, figsize=(14, 3), sharex=True, sharey=True, constrained_layout=True)
axes[0].set_title('Scatter Plot')
axes[1].set_title('Gridded Plot - Default (100x100)')
axes[2].set_title('Gridded Plot - A bit coarser (50x50)')
ptm1.plot(ax=axes[0])
ptm1.plot(plot_type='hist', ax=axes[1], colorbar_kwargs=dict(orientation='horizontal'), cmap='turbo')
ptm1.plot(plot_type='hist', ax=axes[2], colorbar_kwargs=dict(orientation='horizontal'), nx=50, ny=50, cmap='turbo')
plt.show()
As above, but zoomed out a little more and with the mesh turned back on for reference.
fig, axes = plt.subplots(ncols=3, figsize=(10, 5), sharex=True, sharey=True, constrained_layout=True)
axes[0].set_title('Scatter')
axes[1].set_title('Gridded (100x100)')
axes[2].set_title('Gridded (50x50)')
ptm1.plot(ax=axes[0])
ptm1.plot(plot_type='hist', ax=axes[1], colorbar=False, cmap='turbo')
ptm1.plot(plot_type='hist', ax=axes[2], colorbar=False, nx=50, ny=50, cmap='turbo')
for ax in axes:
fv1.plot_mesh(ax=ax, mesh=False)
plt.show()
Plot Variable on Grid¶
We can use the color_by argument to bin up variables and show them on the grid. This is a great way to see things like particle age, or even group or status when the particles are all mixed up
fig, axes = plt.subplots(ncols=3, figsize=(14, 3), sharex=True, sharey=True, constrained_layout=True)
axes[0].set_title('Scatter Plot')
axes[1].set_title('Gridded Plot - Default (100x100)')
axes[2].set_title('Gridded Plot - A bit coarser (50x50)')
ptm1.plot(plot_type='hist', ax=axes[0], colorbar_kwargs=dict(orientation='horizontal'), cmap='turbo', nx=50, ny=50)
ptm1.plot(plot_type='hist', color_by='age', ax=axes[1], colorbar_kwargs=dict(orientation='horizontal'), cmap='turbo', nx=50, ny=50)
ptm1.plot(plot_type='hist', color_by='groupID', ax=axes[2], colorbar_kwargs=dict(orientation='horizontal'), nx=50, ny=50, cmap='turbo')
plt.show()
Plot mass on Grid¶
Some variables like mass and uvw have a third dimension that we need to reduce before we plot them.
print(ptm1['mass'].dims)
('time', 'trajectory', 'mass_consituents')
Let’s say we want to plot the total mass in a grid, across all the mass constituents. In our example case, we only have 1 mass constituent but you may sometimes have a more detailed model that contains several constituents here.
ptm1['total_mass'] = ptm1['mass'].sum(dim='mass_consituents')
Now when we plot, we also need to think about what we want to show. If we just plot the default grid, we’ll end up showing the mean mass of the particles in the grid cell; i.e., what is the mean of the mass that each particle has. But most of the time we actually want the total mass in the cell, so we need to change the statistic from mean (default) to sum.
fig, axes = plt.subplots(ncols=3, figsize=(14, 3), sharex=True, sharey=True, constrained_layout=True)
axes[0].set_title('Density')
axes[1].set_title('Mean of Total Mass Per Particle')
axes[2].set_title('Sum Of Total Mass')
ptm1.plot(plot_type='hist', ax=axes[0], colorbar_kwargs=dict(orientation='horizontal'), cmap='turbo', nx=50, ny=50)
ptm1.plot(plot_type='hist', color_by='total_mass', ax=axes[1], colorbar_kwargs=dict(orientation='horizontal'), cmap='turbo', nx=50, ny=50)
ptm1.plot(plot_type='hist', color_by='total_mass', statistic='sum', ax=axes[2], colorbar_kwargs=dict(orientation='horizontal'), cmap='turbo', nx=50, ny=50)
plt.show()
Plot Interactive¶
We often want to visualise our results interactively to see if the model is behaving as we expect it to.
Similar to hydrodynamic results, we can use plot_interactive to interactively scroll through out timeseries. We need to first enable widget mode with %matplotlib widget, and then we can call the function afterwards.
Any of the options of .plot are available in the interactive mode.
Note: The example below is a static example from another model showing how the interactive mode could look like for a large result
%matplotlib widget
ptm1.plot_interactive(plot_type='hist')
plt.close('all')
%matplotlib inline
This concludes the examples on particle plotting.