🎓 Working with Parcels output#
This tutorial covers the format of the trajectory output exported by Parcels. Parcels does not include advanced analysis or plotting functionality, which users are suggested to write themselves to suit their research goals. Here we provide some starting points to explore the Parcels output files yourself.
For more advanced reading and tutorials on the analysis of Lagrangian trajectories, we recommend checking out the Lagrangian Diagnostics Analysis Cookbook and the project in general.
from datetime import datetime
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
import pyarrow.parquet as pq
import parcels
import parcels.tutorial
First we need to create some Parcels output to analyse. We simulate a set of particles using the setup described in the Delay start tutorial. We will also add some user defined metadata to the output file.
# Load the CopernicusMarine data in the Agulhas region from the example_datasets
ds_fields = parcels.tutorial.open_dataset(
"CopernicusMarine_data_for_Argo_tutorial/data"
).isel(depth=0) # select the first depth level (surface)
# Convert to SGRID-compliant dataset and create FieldSet
fields = {"U": ds_fields["uo"], "V": ds_fields["vo"]}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
# Convert the FieldSet to windowed arrays for better performance
fieldset = fieldset.to_windowed_arrays()
INFO: No depth dimension found in dataset. Assuming surface data and adding a singleton depth dimension.
# Particle locations and initial time
npart = 10 # number of particles to be released
lon = 32 * np.ones(npart)
lat = np.linspace(-32.5, -30.5, npart, dtype=np.float32)
time = ds_fields.time.values[0] + np.arange(0, npart) * np.timedelta64(2, "h")
pset = parcels.ParticleSet(
fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, t=time
)
output_file = parcels.ParticleFile("output.parquet", outputdt=np.timedelta64(2, "h"))
Parcels saves some metadata in the output file with every simulation (Parcels version, CF convention information, etc.). This metadata is just a dictionary which is propogated to the parquet metadata. We are free to manipulate this dictionary to add any custom metadata relevant to our simulation. Here we add a custom metadata field date_created to the output file.
output_file.metadata["date_created"] = datetime.now().isoformat()
output_file.metadata
{'date_created': '2026-10-23T14:05:37.711825'}
To write the metadata to the output_file, we need to add it before running pset.execute() which writes the particleset including the metadata to the output_file.
pset.execute(
parcels.kernels.AdvectionRK2,
runtime=np.timedelta64(48, "h"),
dt=np.timedelta64(5, "m"),
output_file=output_file,
)
Reading the output file#
Note
As of Parcels v4, the default output format is parquet (instead of zarr). The parquet output format is a tabular format, in which every row corresponds to an observation of a particle trajectory. The zarr output format is a multidimensional array format, in which the data is stored in a 2D array with dimensions traj and obs. The parquet format is more compact and faster to read.
However, the parquet format does not support the CF-convention for trajectories data implemented with the NCEI trajectory template. We have implemented a parcels.read_particlefile() function to facilitate reading parquet output files, see more information below.
There are also utilities to convert to v3 ParticleFile output - see later in this tutorial.
Parcels exports output trajectories in parquet format. Files in parquet are stored in tabular data, so each row corresponds to a particle at a given time step, and columns correspond to particle attributes (lon, lat, time, etc.).
The files can be analysed with a wide range of tools, including pandas, polars or the lt toolbox. The latter is specifically designed for the analysis of Lagrangian trajectories, and can be used to compute a wide range of Lagrangian diagnostics- but is still in alpha stage of development.
In Polars, these files can be opened with polars.read_parquet:
df_polars = pl.read_parquet("output.parquet")
print(df_polars)
shape: (205, 5)
┌──────────┬─────┬────────────┬───────────┬─────────────┐
│ t ┆ z ┆ y ┆ x ┆ particle_id │
│ --- ┆ --- ┆ --- ┆ --- ┆ --- │
│ f64 ┆ f32 ┆ f32 ┆ f32 ┆ i64 │
╞══════════╪═════╪════════════╪═══════════╪═════════════╡
│ 0.0 ┆ 0.0 ┆ -32.5 ┆ 32.0 ┆ 0 │
│ 7200.0 ┆ 0.0 ┆ -32.485836 ┆ 31.999071 ┆ 0 │
│ 7200.0 ┆ 0.0 ┆ -32.277779 ┆ 32.0 ┆ 1 │
│ 14400.0 ┆ 0.0 ┆ -32.471802 ┆ 31.997784 ┆ 0 │
│ 14400.0 ┆ 0.0 ┆ -32.262951 ┆ 31.991676 ┆ 1 │
│ … ┆ … ┆ … ┆ … ┆ … │
│ 172800.0 ┆ 0.0 ┆ -31.569975 ┆ 31.236179 ┆ 5 │
│ 172800.0 ┆ 0.0 ┆ -31.509027 ┆ 31.222284 ┆ 6 │
│ 172800.0 ┆ 0.0 ┆ -31.30382 ┆ 31.385489 ┆ 7 │
│ 172800.0 ┆ 0.0 ┆ -31.054222 ┆ 31.595919 ┆ 8 │
│ 172800.0 ┆ 0.0 ┆ -30.785988 ┆ 31.690128 ┆ 9 │
└──────────┴─────┴────────────┴───────────┴─────────────┘
To see our custom metadata field date_created, we need to use pyarrow.parquet
schema = pq.read_schema("output.parquet")
for k, v in schema.metadata.items():
print(k.decode(), "->", v.decode())
date_created -> 2026-10-23T14:05:37.711825
feature_type -> trajectory
Conventions -> CF-1.6/CF-1.7
ncei_template_version -> NCEI_NetCDF_Trajectory_Template_v2.0
parcels_version -> 4.0.3.dev30
parcels_grid_mesh -> SphericalMesh(radius=6366707.019493707)
parcels_kernels -> AdvectionRK2
As you may have noticed above, the t is shown as a float64 (in seconds) in df_polars. That is because polars.read_parquet does not automatically convert the cftime. To handle this, we also provide a helper function parcels.read_particlefile() to read ParticleFiles, which does automatically convert the cftime.
df_particles = parcels.read_particlefile("output.parquet")
print(df_particles)
shape: (205, 5)
┌─────────────────────┬─────┬────────────┬───────────┬─────────────┐
│ t ┆ z ┆ y ┆ x ┆ particle_id │
│ --- ┆ --- ┆ --- ┆ --- ┆ --- │
│ datetime[ns] ┆ f32 ┆ f32 ┆ f32 ┆ i64 │
╞═════════════════════╪═════╪════════════╪═══════════╪═════════════╡
│ 2024-01-01 00:00:00 ┆ 0.0 ┆ -32.5 ┆ 32.0 ┆ 0 │
│ 2024-01-01 02:00:00 ┆ 0.0 ┆ -32.485836 ┆ 31.999071 ┆ 0 │
│ 2024-01-01 02:00:00 ┆ 0.0 ┆ -32.277779 ┆ 32.0 ┆ 1 │
│ 2024-01-01 04:00:00 ┆ 0.0 ┆ -32.471802 ┆ 31.997784 ┆ 0 │
│ 2024-01-01 04:00:00 ┆ 0.0 ┆ -32.262951 ┆ 31.991676 ┆ 1 │
│ … ┆ … ┆ … ┆ … ┆ … │
│ 2024-01-03 00:00:00 ┆ 0.0 ┆ -31.569975 ┆ 31.236179 ┆ 5 │
│ 2024-01-03 00:00:00 ┆ 0.0 ┆ -31.509027 ┆ 31.222284 ┆ 6 │
│ 2024-01-03 00:00:00 ┆ 0.0 ┆ -31.30382 ┆ 31.385489 ┆ 7 │
│ 2024-01-03 00:00:00 ┆ 0.0 ┆ -31.054222 ┆ 31.595919 ┆ 8 │
│ 2024-01-03 00:00:00 ┆ 0.0 ┆ -30.785988 ┆ 31.690128 ┆ 9 │
└─────────────────────┴─────┴────────────┴───────────┴─────────────┘
Trajectory data structure#
The output dataset used here contains 10 particles and 25 observations. Not every particle has 25 observations however; since we released particles at different times some particle trajectories are shorter than others.
for traj in df_particles.partition_by("particle_id", maintain_order=True):
time_origin = pd.Timestamp(fieldset.time_interval.left).to_pydatetime()
time_in_hour = (traj["t"] - time_origin).dt.total_hours()
traj_id = traj["particle_id"][0]
print(f"Particle {traj_id}: " + "".join(f"{int(t):2d} " for t in time_in_hour))
Particle 0: 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 1: 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 2: 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 3: 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 4: 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 5: 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 6: 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 7: 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 8: 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Particle 9: 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
Note how the first observation occurs at a different time for each trajectory.
Plotting trajectories#
Parcels output consists of particle trajectories through time and space. An important way to explore patterns in this information is to draw the trajectories in space. The trajan package can be used to quickly plot Parcels results, but users are encouraged to create their own figures, for example by using the comprehensive matplotlib library. Here we show a basic setup on how to process the Parcels output into trajectory plots and animations.
Warning
Trajan is not yet compatible with the parquet output format, but we are working on a solution to this.
Since the parquet output format is tabular, a simple plot of the longitude and latitude of all particles will show one continuous line:
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(df_particles["x"], df_particles["y"], ".-")
plt.show()
To show the trajectories of individual particles, we must group the data by trajectory and plot each trajectory separately:
fig, ax = plt.subplots(figsize=(5, 3))
for traj in df_particles.partition_by("particle_id", maintain_order=True):
traj_id = traj["particle_id"][0]
ax.plot(traj["x"], traj["y"], ".-", label=f"P{traj_id}")
ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5), borderaxespad=0.0)
plt.tight_layout()
plt.show()
However, if you want to plot only the locations at a certain time step, you can simply select the data at that time step:
time_step = np.timedelta64(18, "h")
time_to_plot = fieldset.time_interval.left + time_step
particles = df_particles.filter(pl.col("t") == pl.lit(time_to_plot))
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(particles["x"], particles["y"], "o")
title_time = pd.to_datetime(time_to_plot).strftime("%Y-%m-%d %H:%M:%S")
ax.set_title(f"Particle locations at {title_time}")
plt.show()
Or, if you want to plot the particles a certain amount of time after they were released, you can select the data based on the time since release:
time_step = np.timedelta64(18, "h")
particles = df_particles.filter(
(pl.col("t") - pl.col("t").min().over("particle_id")) == pl.lit(time_step)
)
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(particles["x"], particles["y"], "o")
ax.set_title(f"Particle locations {time_step} after their release")
plt.show()
We can plot the distance travelled as a function of absolute and relative time:
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
for traj in df_particles.partition_by("particle_id", maintain_order=True):
distance = np.sqrt(
(traj["x"] - traj["x"][0]) ** 2 + (traj["y"] - traj["y"][0]) ** 2
)
ax[0].plot(traj["t"], distance, ".-", label=f"P{traj['particle_id'][0]}")
rel_time = (traj["t"] - traj["t"][0]).dt.total_hours()
ax[1].plot(rel_time, distance, ".-", label=f"P{traj['particle_id'][0]}")
ax[0].set_xlabel("Date")
ax[0].set_ylabel("Distance travelled [degrees]")
ax[0].tick_params(axis="x", labelrotation=45)
ax[1].set_xlabel("Time since release [hours]")
ax[1].set_ylabel("Distance travelled [degrees]")
plt.tight_layout()
plt.show()
Animations#
Trajectory plots like the ones above can become very cluttered for large sets of particles. To better see patterns, it’s a good idea to create an animation in time and space. To do this, matplotlib offers an animation package. Here we show how to use the FuncAnimation class to animate Parcels trajectory data.
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib
from matplotlib.animation import FuncAnimation
from matplotlib.collections import LineCollection
# for interactive display of animation
plt.rcParams["animation.html"] = "jshtml"
time_step = np.timedelta64(2, "h") # time step for animation frames
# Build the time axis for the animation.
timerange = np.arange(
np.nanmin(df_particles["t"]),
np.nanmax(df_particles["t"]) + time_step,
time_step,
)
# figure setup
fig, ax = plt.subplots(figsize=(6, 5), subplot_kw={"projection": ccrs.PlateCarree()})
ax.set_xlim(30, 33)
ax.set_xticks(np.arange(30, 33.5, 0.5))
ax.set_xlabel("Longitude (deg E)")
ax.set_ylim(-33, -30)
ax.set_yticks(ticks=np.arange(-33, -29.5, 0.5))
ax.set_yticklabels(np.arange(33, 29.5, -0.5).astype(str))
ax.set_ylabel("Latitude (deg S)")
ax.coastlines()
ax.add_feature(cfeature.LAND)
cmap = matplotlib.colormaps["tab20b"]
trails = LineCollection([], linewidths=0.6, alpha=0.3)
ax.add_collection(trails)
particle_ids = df_particles["particle_id"].unique().to_list()
# Convert the trajectory data to a time-by-particle array and assign colors.
colors = np.asarray(
[cmap(i / max(len(particle_ids) - 1, 1)) for i in range(len(particle_ids))]
)
# Create arrays to hold the x and y positions of each particle at each time step.
x = np.full((len(particle_ids), len(timerange)), np.nan)
y = np.full((len(particle_ids), len(timerange)), np.nan)
traj = df_particles.with_columns(
pl.col("particle_id")
.replace(particle_ids, range(len(particle_ids)))
.alias("p_idx"),
pl.int_range(pl.len()).alias("row_idx"),
)
# Map each observation to its nearest animation time index.
traj = traj.with_columns(
pl.Series("t_idx", np.searchsorted(timerange, traj["t"].to_numpy())).alias("t_idx")
)
x[traj["p_idx"].to_numpy(), traj["t_idx"].to_numpy()] = traj["x"].to_numpy()
y[traj["p_idx"].to_numpy(), traj["t_idx"].to_numpy()] = traj["y"].to_numpy()
# Plot first timestep
scatter = ax.scatter(x[:, 0], y[:, 0], s=10, c=colors, zorder=2)
# Set initial title
t_str = pd.to_datetime(timerange[0]).strftime("%Y-%m-%d %H:%M:%S")
title = ax.set_title(f"Particles at t = {t_str}")
def animate(i):
t_str = pd.to_datetime(timerange[i]).strftime("%Y-%m-%d %H:%M:%S")
title.set_text(f"Particles at t = {t_str}")
scatter.set_offsets(np.column_stack((x[:, i], y[:, i])))
trail_length = min(10, i) # trails have max length of 10 time steps
start = max(0, i - trail_length)
x_slice = x[:, start : i + 1]
y_slice = y[:, start : i + 1]
trails.set_segments(np.dstack((x_slice, y_slice)))
trails.set_color(colors)
# Create animation
anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)
plt.close(fig)
anim
/home/docs/checkouts/readthedocs.org/user_builds/parcels/checkouts/latest/.pixi/envs/docs/lib/python3.15/site-packages/cartopy/io/__init__.py:263: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_land.zip
warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/docs/checkouts/readthedocs.org/user_builds/parcels/checkouts/latest/.pixi/envs/docs/lib/python3.15/site-packages/cartopy/io/__init__.py:263: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
warnings.warn(f'Downloading: {url}', DownloadWarning)
Using v3 output#
We also have some tooling for converting v4 trajectory output to the v3 Zarr output. This is meant as a way to reduce friction for those not wanting to update their plotting scripts, however we recommend writing plotting scripts that work with the new v4 output.
At the moment, this tooling is not lazy (i.e., it eagerly loads all data into memory in order to do the transformation) so your mileage may vary if you try to convert larger datasets. If you would like to contribute to the performance of this utility, please open an issue or PR.
import xarray as xr
parcels.particlefile_to_v3_zarr("output.parquet", "output.zarr")
ds = xr.open_dataset("output.zarr")
ds
<xarray.Dataset> Size: 5kB
Dimensions: (trajectory: 10, obs: 25)
Coordinates:
* trajectory (trajectory) int64 80B 0 1 2 3 4 5 6 7 8 9
* obs (obs) int64 200B 0 1 2 3 4 5 6 7 8 ... 17 18 19 20 21 22 23 24
Data variables:
lat (trajectory, obs) float32 1kB ...
lon (trajectory, obs) float32 1kB ...
time (trajectory, obs) datetime64[ns] 2kB ...
z (trajectory, obs) float32 1kB ...
Attributes:
date_created: 2026-10-23T14:05:37.711825
feature_type: trajectory
Conventions: CF-1.6/CF-1.7
ncei_template_version: NCEI_NetCDF_Trajectory_Template_v2.0
parcels_version: 4.0.3.dev30
parcels_grid_mesh: SphericalMesh(radius=6366707.019493707)
parcels_kernels: AdvectionRK2