Skip to content

Commit f644537

Browse files
erikvansebilleerikvansebille
authored
Merge pull request #2479 from Parcels-code/tutorial_summingfields
Adding a tutorial on how to manipulate fields using xarray operations
2 parents 3025125 + 2801d59 commit f644537

File tree

3 files changed

+193
-0
lines changed

3 files changed

+193
-0
lines changed

‎docs/user_guide/examples/explanation_kernelloop.md‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,10 @@ Besides having commutable Kernels, the main advantage of this implementation is
4242

4343
Below is a simple example of some particles at the surface of the ocean. We create an idealised zonal wind flow that will "push" a particle that is already affected by the surface currents. The Kernel loop ensures that these two forces act at the same time and location.
4444

45+
```{note}
46+
Advecting particles by a combination of flow fields can be done with two separate kernels, as is shown below. However, it can also be done by summing the fields directly using `xarray` operations - provided the fields are defined on the same grid and have compatible dimensions. See the [manipulating field data tutorial](tutorial_manipulating_field_data.ipynb) for an example of that approach.
47+
```
48+
4549
```{code-cell}
4650
:tags: [hide-output]
4751
import matplotlib.pyplot as plt

‎docs/user_guide/examples/tutorial_manipulating_field_data.ipynb‎

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "0",
6+
"metadata": {},
7+
"source": [
8+
"# 🖥️ Manipulating Field data"
9+
]
10+
},
11+
{
12+
"cell_type": "markdown",
13+
"id": "1",
14+
"metadata": {},
15+
"source": [
16+
"Because Parcels `Field` objects are built on top of `xarray` DataArrays, 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 examples of how to leverage that power."
17+
]
18+
},
19+
{
20+
"cell_type": "markdown",
21+
"id": "2",
22+
"metadata": {},
23+
"source": [
24+
"## Summing Fields\n",
25+
"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.\n",
26+
"\n",
27+
"We base this tutorial on the example provided in the [Parcels Kernel loop explanation](explanation_kernelloop.md). There, we combined two kernels to simulate the joint effect of currents and winds.\n",
28+
"\n",
29+
"However, we can also do that in one kernel, by combining the fields directly, using `xarray` operations.\n",
30+
"\n",
31+
"We start with loading the necessary libraries and data"
32+
]
33+
},
34+
{
35+
"cell_type": "code",
36+
"execution_count": null,
37+
"id": "3",
38+
"metadata": {},
39+
"outputs": [],
40+
"source": [
41+
"import matplotlib.pyplot as plt\n",
42+
"import numpy as np\n",
43+
"import xarray as xr\n",
44+
"\n",
45+
"import parcels\n",
46+
"\n",
47+
"# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n",
48+
"example_dataset_folder = parcels.download_example_dataset(\n",
49+
" \"CopernicusMarine_data_for_Argo_tutorial\"\n",
50+
")\n",
51+
"\n",
52+
"ds_fields = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n",
53+
"ds_fields.load() # load the dataset into memory\n",
54+
"\n",
55+
"# Create an idealised wind field and add it to the dataset\n",
56+
"tdim, ydim, xdim = (\n",
57+
" len(ds_fields.time),\n",
58+
" len(ds_fields.latitude),\n",
59+
" len(ds_fields.longitude),\n",
60+
")\n",
61+
"ds_fields[\"UWind\"] = xr.DataArray(\n",
62+
" data=np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values)[None, :, None],\n",
63+
" coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n",
64+
")\n",
65+
"\n",
66+
"ds_fields[\"VWind\"] = xr.DataArray(\n",
67+
" data=np.zeros((tdim, ydim, xdim)),\n",
68+
" coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n",
69+
")"
70+
]
71+
},
72+
{
73+
"cell_type": "markdown",
74+
"id": "4",
75+
"metadata": {},
76+
"source": [
77+
"Now here comes the trick: we can simply sum the fields together using `xarray` operations"
78+
]
79+
},
80+
{
81+
"cell_type": "code",
82+
"execution_count": null,
83+
"id": "5",
84+
"metadata": {},
85+
"outputs": [],
86+
"source": [
87+
"# Combine ocean currents and wind fields\n",
88+
"ds_fields[\"U\"] = ds_fields[\"uo\"] + ds_fields[\"UWind\"]\n",
89+
"ds_fields[\"V\"] = ds_fields[\"vo\"] + ds_fields[\"VWind\"]"
90+
]
91+
},
92+
{
93+
"cell_type": "markdown",
94+
"id": "6",
95+
"metadata": {},
96+
"source": [
97+
"```{note}\n",
98+
"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.\n",
99+
"```"
100+
]
101+
},
102+
{
103+
"cell_type": "markdown",
104+
"id": "7",
105+
"metadata": {},
106+
"source": [
107+
"We can then run the same simulation as in the [Kernel loop explanation tutorial](explanation_kernelloop.md), 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."
108+
]
109+
},
110+
{
111+
"cell_type": "code",
112+
"execution_count": null,
113+
"id": "8",
114+
"metadata": {
115+
"tags": [
116+
"hide-output"
117+
]
118+
},
119+
"outputs": [],
120+
"source": [
121+
"fields = {\n",
122+
" \"U\": ds_fields[\"U\"],\n",
123+
" \"V\": ds_fields[\"V\"],\n",
124+
"}\n",
125+
"ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)\n",
126+
"fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)\n",
127+
"\n",
128+
"npart = 10\n",
129+
"z = np.repeat(ds_fields.depth[0].values, npart)\n",
130+
"lons = np.repeat(31, npart)\n",
131+
"lats = np.linspace(-32.5, -30.5, npart)\n",
132+
"\n",
133+
"pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, lat=lats, lon=lons)\n",
134+
"output_file = parcels.ParticleFile(\n",
135+
" store=\"summed_advection_wind.zarr\", outputdt=np.timedelta64(6, \"h\")\n",
136+
")\n",
137+
"pset.execute(\n",
138+
" [parcels.kernels.AdvectionRK4],\n",
139+
" runtime=np.timedelta64(5, \"D\"),\n",
140+
" dt=np.timedelta64(1, \"h\"),\n",
141+
" output_file=output_file,\n",
142+
")"
143+
]
144+
},
145+
{
146+
"cell_type": "markdown",
147+
"id": "9",
148+
"metadata": {},
149+
"source": [
150+
"We can then plot the trajectories, and confirm that they particles move the same as in the [Kernel loop explanation tutorial](explanation_kernelloop.md), where we combined the effects of currents and winds using two separate kernels."
151+
]
152+
},
153+
{
154+
"cell_type": "code",
155+
"execution_count": null,
156+
"id": "10",
157+
"metadata": {},
158+
"outputs": [],
159+
"source": [
160+
"# Plot the resulting particle trajectories overlapped for both cases\n",
161+
"summed_advection_wind = xr.open_zarr(\"summed_advection_wind.zarr\")\n",
162+
"plt.plot(summed_advection_wind.lon.T, summed_advection_wind.lat.T, \"-\")\n",
163+
"plt.show()"
164+
]
165+
}
166+
],
167+
"metadata": {
168+
"kernelspec": {
169+
"display_name": "docs",
170+
"language": "python",
171+
"name": "python3"
172+
},
173+
"language_info": {
174+
"codemirror_mode": {
175+
"name": "ipython",
176+
"version": 3
177+
},
178+
"file_extension": ".py",
179+
"mimetype": "text/x-python",
180+
"name": "python",
181+
"nbconvert_exporter": "python",
182+
"pygments_lexer": "ipython3",
183+
"version": "3.14.2"
184+
}
185+
},
186+
"nbformat": 4,
187+
"nbformat_minor": 5
188+
}

‎docs/user_guide/index.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ examples/tutorial_nemo_curvilinear.ipynb
2929
examples/tutorial_nemo_3D.ipynb
3030
examples/tutorial_velocityconversion.ipynb
3131
examples/tutorial_nestedgrids.ipynb
32+
examples/tutorial_summingfields.ipynb
3233
```
3334

3435
<!-- examples/documentation_indexing.ipynb -->

0 commit comments

Comments
 (0)