🖥️ Manipulating Field data#
Because parcels.FieldSet objects are built on top of xarray DataSets, you can leverage the powerful data manipulation capabilities of xarray to preprocess or modify your field data before using it in particle simulations. This tutorial provides some common examples of how to leverage that power - however we highly recommend also reading the Xarray documentation!
Let’s start with some imports.
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import parcels
import parcels.tutorial
Looping field data and overriding the time dimension#
Sometimes we want to loop field data (e.g., use one year of data again ten times), or we want to completely override the time dimension.
In previous versions of Parcels (pre version 4) this was done via Parcels function parameters timestamps and time_periodic. In this version of Parcels, we simply just manipulate the Xarray objects directly before passing the objects to Parcels.
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.2Here we’re dealing with daily data for from 2024-01-01 to 2024-02-01 (32 observations in total).
Let’s loop this data for a total of 90 days.
N_TIME = ds_fields.time.size
N_TIME_TARGET = 90
time_indices = np.arange(N_TIME_TARGET) % N_TIME
new_time_coordinate = np.arange(
np.datetime64("2000-01-01"),
np.datetime64("2000-01-01") + N_TIME_TARGET * np.timedelta64(1, "D"),
np.timedelta64(1, "D"),
).astype("datetime64[ns]")
print(time_indices)
[ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
24 25 26 27 28 29 30 31 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 0 1 2 3 4 5 6 7
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25]
ds_periodic = ds_fields.isel(time=time_indices)
ds_periodic["time"] = new_time_coordinate
# Now we have our periodic Dataset, let's take a look at the U velocity at a given position
# to verify
ds_periodic["uo"].isel(depth=0, latitude=0, longitude=0).plot()
[<matplotlib.lines.Line2D at 0x74628b8f70e0>]
From this example, its easy to see how we can completely override the time dimension. For example, by adding a linear offset to the coordinates to make the data for February (and the first few days of March).
new_time_coordinate = ds_fields.time.values + np.timedelta64(31, "D")
ds_offset = ds_fields.copy()
ds_offset["time"] = new_time_coordinate
ds_offset
<xarray.Dataset> Size: 19MB
Dimensions: (time: 32, depth: 41, latitude: 37, longitude: 25)
Coordinates:
* time (time) datetime64[ns] 256B 2024-02-01 2024-02-02 ... 2024-03-03
* 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.2Summing Fields#
In some applications, you may want to sum multiple fields together to create a combined effect on particle movement. For example, you might want to combine ocean currents with wind-driven surface drift. This tutorial demonstrates how to sum multiple fields in Parcels and visualize the resulting particle trajectories.
We base this tutorial on the example provided in the Parcels Kernel loop explanation. There, we combined two kernels to simulate the joint effect of currents and winds.
However, we can also do that in one kernel, by combining the fields directly, using xarray operations.
We start with loading the necessary libraries and data
# Load the CopernicusMarine data in the Agulhas region from the example_datasets
ds_fields = parcels.tutorial.open_dataset(
"CopernicusMarine_data_for_Argo_tutorial/data"
)
# Create an idealised wind field and add it to the dataset
ydim, xdim = len(ds_fields.latitude), len(ds_fields.longitude)
ds_fields["UWind"] = xr.DataArray(
data=0.5
* np.ones((ydim, xdim))
* np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[:, None],
coords=[ds_fields.latitude, ds_fields.longitude],
)
ds_fields["VWind"] = xr.DataArray(
data=np.zeros((ydim, xdim)),
coords=[ds_fields.latitude, ds_fields.longitude],
)
Now here comes the trick: we can simply sum the fields together using xarray operations
# Combine ocean currents and wind fields
ds_fields["U"] = ds_fields["uo"] + ds_fields["UWind"]
ds_fields["V"] = ds_fields["vo"] + ds_fields["VWind"]
Note
Combining fields in this way assumes that the fields are defined on the same grid and have compatible dimensions. Ensure that the fields you are summing are aligned correctly to avoid unexpected results.
We can then run the same simulation as in the Kernel loop explanation tutorial, but now using the combined fields. We only need the default advection kernel now, since the effects of both currents and winds are already included in the summed fields.
fields = {
"U": ds_fields["U"],
"V": ds_fields["V"],
}
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()
npart = 10
lons = np.repeat(32.2, npart)
lats = np.linspace(-32.5, -30.5, npart)
pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, y=lats, x=lons)
output_file = parcels.ParticleFile(
path="summed_advection_wind.parquet",
outputdt=np.timedelta64(6, "h"),
)
pset.execute(
[parcels.kernels.AdvectionRK2],
runtime=np.timedelta64(5, "D"),
dt=np.timedelta64(1, "h"),
output_file=output_file,
)
We can then plot the trajectories, and confirm that they particles move the same as in the Kernel loop explanation tutorial, where we combined the effects of currents and winds using two separate kernels.
# Plot the resulting particle trajectories overlapped for both cases
summed_advection_wind = parcels.read_particlefile("summed_advection_wind.parquet")
fig, ax = plt.subplots(figsize=(5, 3))
for traj in summed_advection_wind.partition_by("particle_id", maintain_order=True):
ax.plot(traj["x"], traj["y"], "-")
plt.show()
