🎓 Quickstart tutorial#
Welcome to the Parcels quickstart tutorial, in which we will go through all the necessary steps to run a simulation. The code in this notebook can be used as a starting point to run Parcels in your own environment. Along the way we will familiarize ourselves with some specific classes and methods. If you are ever confused about one of these and want to read more, we have a concepts overview discussing them in more detail. Let’s dive in!
Imports#
Parcels depends on xarray, expecting inputs in the form of xarray.Dataset. Output files can be read with Polars.
import numpy as np
import xarray as xr
import polars as pl
import parcels
import parcels.tutorial
Input flow fields: FieldSet#
A Parcels simulation of Lagrangian trajectories of virtual particles requires two inputs; the first is a set of hydrodynamics fields in which the particles are tracked. Here we provide an example using a subset of the Global Ocean Physics Reanalysis from the Copernicus Marine Service.
ds_fields = parcels.tutorial.open_dataset("CopernicusMarine_data_for_Argo_tutorial/data")
ds_fields
<xarray.Dataset> Size: 19MB
Dimensions: (time: 32, depth: 41, latitude: 37, longitude: 25)
Coordinates:
* time (time) datetime64[ns] 256B 2024-01-01 2024-01-02 ... 2024-02-01
* depth (depth) float32 164B 0.494 1.541 2.646 ... 1.942e+03 2.225e+03
* latitude (latitude) float32 148B -33.0 -32.92 -32.83 ... -30.08 -30.0
* longitude (longitude) float32 100B 31.0 31.08 31.17 ... 32.83 32.92 33.0
Data variables:
uo (time, depth, latitude, longitude) float32 5MB dask.array<chunksize=(32, 41, 37, 25), meta=np.ndarray>
vo (time, depth, latitude, longitude) float32 5MB dask.array<chunksize=(32, 41, 37, 25), meta=np.ndarray>
so (time, depth, latitude, longitude) float32 5MB dask.array<chunksize=(32, 41, 37, 25), meta=np.ndarray>
thetao (time, depth, latitude, longitude) float32 5MB dask.array<chunksize=(32, 41, 37, 25), meta=np.ndarray>
Attributes:
Conventions: CF-1.8
credit: E.U. Copernicus Marine Service Information (CM...
institution: Mercator Ocean International
title: daily mean fields from Global Ocean Physics An...
source: MOI GLO12
contact: https://marine.copernicus.eu/contact
references: http://marine.copernicus.eu
producer: CMEMS - Global Monitoring and Forecasting Centre
copernicusmarine_version: 2.2.2As we can see, the reanalysis dataset contains eastward velocity uo, northward velocity vo, potential temperature
(thetao) and salinity (so) fields.
These hydrodynamic fields need to be stored in a parcels.FieldSet object. Parcels provides tooling to parse many types
of models or observations into such a parcels.FieldSet object. This is done in a two-step approach.
First, we convert the dataset into an SGRID-compliant dataset, for example by using a version of parcels.convert.<MODEL>_to_sgrid(). Then, we create the parcels.FieldSet from the SGRID-compliant dataset using parcels.FieldSet.from_sgrid_conventions().
Below, we use a combination of parcels.convert.copernicusmarine_to_sgrid() and parcels.FieldSet.from_sgrid_conventions(), providing the names of the velocity fields in the dataset in the dictionary fields:
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)
Now, in order to improve performance, we can convert the parcels.FieldSet to windowed arrays. This is especially useful for large datasets with many timeslices, as it allows Parcels to load only the necessary timeslices into memory during the simulation. For more information about squeezing performance out of Parcels, see the performance tutorial.
# Convert the FieldSet to windowed arrays for better performance
fieldset = fieldset.to_windowed_arrays()
You can inspect the parcels.FieldSet object with the describe method in order to see which parcels.Fields are included, and which grid and interpolation method is used for each field. This also gives information on the type of mesh and the time interval of the parcels.FieldSet:
fieldset.describe()
| Name | Type | Grid number | Interp method / value | Parcels backend |
|:-------|:------------|--------------:|:------------------------|:------------------|
| U | Field | 0 | XLinear(...) | WindowedArray |
| V | Field | 0 | XLinear(...) | WindowedArray |
| UV | VectorField | 0 | XLinear_Velocity(...) | - |
mesh: SphericalMesh(radius=6366707.019493707)
time interval: (np.datetime64('2024-01-01T00:00:00.000000000'), np.datetime64('2024-02-01T00:00:00.000000000'))
The subset contains a region of the Agulhas current along the southeastern coast of Africa:
temperature = ds_fields.isel(time=0, depth=0).thetao.plot(cmap="magma")
velocity = ds_fields.isel(time=0, depth=0).plot.quiver(x="longitude", y="latitude", u="uo", v="vo")
Input virtual particles: ParticleSet#
Now that we have created a parcels.FieldSet object from the hydrodynamic data, we need to provide our second input:
the virtual particles for which we will calculate the trajectories.
We need to create a parcels.ParticleSet object with the particles’ initial time and position. The parcels.ParticleSet
object also needs to know about the parcels.FieldSet in which the particles “live”. Finally, we need to specify the type of
parcels.ParticleClass we want to use. The default particles have t, z, y, and x, but you can easily add
other parcels.Variables such as size, temperature, or age to create your own particles to mimic plastic or an ARGO float.
# Particle locations and initial time
npart = 10 # number of particles to be released
# release particles in a line along a meridian
lat = np.linspace(-32.5, -30.5, npart)
lon = np.repeat(32, npart)
time = np.repeat(ds_fields.time.values[0], npart) # at initial time of input data
z = np.repeat(0.5, npart) # at 0.5 meter depth (first level of the input data)
pset = parcels.ParticleSet(
fieldset=fieldset, pclass=parcels.Particle, t=time, z=z, y=lat, x=lon
)
And you can plot the particles on top of the temperature and velocity field:
temperature = ds_fields.isel(time=0, depth=0).thetao.plot(cmap="magma")
velocity = ds_fields.isel(time=0, depth=0).plot.quiver(x="longitude", y="latitude", u="uo", v="vo")
ax = temperature.axes
ax.scatter(lon, lat, s=40, c='w', edgecolors='r');
If you also want to sample the initial value of a field, you can do so by calling the field with the parcels.ParticleSet as an argument. For example, in the code above you would add one line (assuming that the ParticleSet has a variable temperature) to sample the initial temperature values, before calling pset.execute(...).
pset.temperature = fieldset.thetao[pset]
See the sampling tutorial for more information on initial sampling of fields.
Compute: Kernel#
After setting up the input data and particle start locations and times, we need to specify what calculations to do with
the particles. These calculations, or numerical integrations, will be performed by what we call a kernel function, operating on
all particles in the ParticleSet. The most common calculation is the advection of particles through the velocity field.
Parcels comes with a number of common parcels.kernels, from which we will use the Runge-Kutta advection kernel parcels.kernels.AdvectionRK2:
kernels = [parcels.kernels.AdvectionRK2]
Prepare output: ParticleFile#
Before starting the simulation, we must define where and how frequent we want to write the output of our simulation.
We can define this in a parcels.ParticleFile object:
output_file = parcels.ParticleFile("output-quickstart.parquet", outputdt=np.timedelta64(1, "h"))
The output files are in .parquet format, which can be read by Polars.
See the Parcels output tutorial for more information on the parquet format. We want to choose
the outputdt argument so that it captures the smallest timescales of our interest.
Run Simulation: ParticleSet.execute()#
Finally, we can run the simulation by executing the ParticleSet using the specified list of kernels. This is done using the parcels.ParticleSet.execute() method.
Additionally, we need to specify:
the
runtime: for how long we want to simulate particles.the
dt: the timestep with which to perform the numerical integration in thekernels. Depending on the numerical integration scheme, the accuracy of our simulation will depend ondt. Read this notebook to learn more about numerical accuracy.
pset.execute(
kernels,
runtime=np.timedelta64(1, "D"),
dt=np.timedelta64(5, "m"),
output_file=output_file,
)
Read output#
To start analysing the trajectories computed by Parcels, we can open the parcels.ParticleFile using the parcels.read_particlefile() utility, which itself uses polars:
df = parcels.read_particlefile("output-quickstart.parquet")
df
| t | z | y | x | particle_id |
|---|---|---|---|---|
| datetime[ns] | f32 | f32 | f32 | i64 |
| 2024-01-01 00:00:00 | 0.5 | -32.5 | 32.0 | 0 |
| 2024-01-01 00:00:00 | 0.5 | -32.277779 | 32.0 | 1 |
| 2024-01-01 00:00:00 | 0.5 | -32.055557 | 32.0 | 2 |
| 2024-01-01 00:00:00 | 0.5 | -31.833334 | 32.0 | 3 |
| 2024-01-01 00:00:00 | 0.5 | -31.611111 | 32.0 | 4 |
| … | … | … | … | … |
| 2024-01-02 00:00:00 | 0.5 | -31.408091 | 31.510376 | 5 |
| 2024-01-02 00:00:00 | 0.5 | -31.368313 | 31.389574 | 6 |
| 2024-01-02 00:00:00 | 0.5 | -31.209724 | 31.438147 | 7 |
| 2024-01-02 00:00:00 | 0.5 | -31.002275 | 31.592552 | 8 |
| 2024-01-02 00:00:00 | 0.5 | -30.758671 | 31.720886 | 9 |
The file contains 250 rows: 25 observations for the 10 particle trajectories. The working with Parcels output tutorial provides more detail about the dataset and how to analyse it.
Let’s verify that Parcels has computed the advection of the virtual particles!
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
# convert time to matplotlib date format for plotting
time_values = mdates.date2num(df["t"].to_list())
# plot positions and colour particles by time
scatter = plt.scatter(df['x'], df['y'], c=time_values)
plt.scatter(df['x'][:npart], df['y'][:npart], facecolors="none", edgecolors='r') # starting positions
plt.scatter(lon, lat, facecolors="none", edgecolors='r') # starting positions
plt.xlim(31,33)
plt.ylabel("Latitude [deg N]")
plt.ylim(-33,-30)
colorbar = plt.colorbar(scatter, label="Time")
colorbar.ax.yaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d %H:%M"))
plt.show()
That looks good! The virtual particles released in a line along the 32nd meridian (dark blue) have been advected by the flow field.
Running a simulation backwards in time#
Now that we know how to run a simulation, we can easily run another and change one of the settings. We can trace back
the particles from their current to their original position by running the simulation backwards in time. To do so, we
can simply make dt < 0.
Note
We have not edited the ParticleSet, which means that the new simulation will start with the particles at their current
location!
# set up output file
output_file = parcels.ParticleFile("output-backwards.parquet", outputdt=np.timedelta64(1, "h"))
# execute simulation in backwards time
pset.execute(
kernels,
runtime=np.timedelta64(1, "D"),
dt=-np.timedelta64(5, "m"),
output_file=output_file,
)
When we check the output, we can see that the particles have returned to their original position!
df_back = parcels.read_particlefile("output-backwards.parquet")
# convert time to matplotlib date format for plotting
time_values = mdates.date2num(df_back["t"].to_list())
scatter = plt.scatter(df_back['x'], df_back['y'], c=time_values)
particles_at_max_time = df_back.filter(pl.col("t") == df_back["t"].max())
plt.scatter(particles_at_max_time['x'], particles_at_max_time['y'], facecolors="none", edgecolors='r') # starting positions
plt.xlabel("Longitude [deg E]")
plt.xlim(31,33)
plt.ylabel("Latitude [deg N]")
plt.ylim(-33,-30)
colorbar = plt.colorbar(scatter, label="Time")
colorbar.ax.yaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d %H:%M"))
plt.show()
Using RK2 advection, the final positions are equal to the original positions with an accuracy of 5 decimals:
# testing that final location == original location
particles_at_min_time = df_back.filter(pl.col("t") == df_back["t"].min())
np.testing.assert_almost_equal(particles_at_min_time["y"], lat, decimal=5)
np.testing.assert_almost_equal(particles_at_min_time['x'], lon, decimal=5)