Adding a tutorial on how to manipulate fields using xarray operations

docs/user_guide/examples/explanation_kernelloop.md

@@ -42,6 +42,10 @@ Besides having commutable Kernels, the main advantage of this implementation is
 
 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.
 
+```{note}
+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.
+```
+
 ```{code-cell}
 :tags: [hide-output]
 import matplotlib.pyplot as plt

docs/user_guide/examples/tutorial_manipulating_field_data.ipynb

@@ -0,0 +1,188 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# 🖥️ Manipulating Field data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2",
+   "metadata": {},
+   "source": [
+    "## Summing Fields\n",
+    "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",
+    "\n",
+    "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",
+    "\n",
+    "However, we can also do that in one kernel, by combining the fields directly, using `xarray` operations.\n",
+    "\n",
+    "We start with loading the necessary libraries and data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "\n",
+    "import parcels\n",
+    "\n",
+    "# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n",
+    "example_dataset_folder = parcels.download_example_dataset(\n",
+    "    \"CopernicusMarine_data_for_Argo_tutorial\"\n",
+    ")\n",
+    "\n",
+    "ds_fields = xr.open_mfdataset(f\"{example_dataset_folder}/*.nc\", combine=\"by_coords\")\n",
+    "ds_fields.load()  # load the dataset into memory\n",
+    "\n",
+    "# Create an idealised wind field and add it to the dataset\n",
+    "tdim, ydim, xdim = (\n",
+    "    len(ds_fields.time),\n",
+    "    len(ds_fields.latitude),\n",
+    "    len(ds_fields.longitude),\n",
+    ")\n",
+    "ds_fields[\"UWind\"] = xr.DataArray(\n",
+    "    data=np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values)[None, :, None],\n",
+    "    coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n",
+    ")\n",
+    "\n",
+    "ds_fields[\"VWind\"] = xr.DataArray(\n",
+    "    data=np.zeros((tdim, ydim, xdim)),\n",
+    "    coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4",
+   "metadata": {},
+   "source": [
+    "Now here comes the trick: we can simply sum the fields together using `xarray` operations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Combine ocean currents and wind fields\n",
+    "ds_fields[\"U\"] = ds_fields[\"uo\"] + ds_fields[\"UWind\"]\n",
+    "ds_fields[\"V\"] = ds_fields[\"vo\"] + ds_fields[\"VWind\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "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",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {
+    "tags": [
+     "hide-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "fields = {\n",
+    "    \"U\": ds_fields[\"U\"],\n",
+    "    \"V\": ds_fields[\"V\"],\n",
+    "}\n",
+    "ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)\n",
+    "fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)\n",
+    "\n",
+    "npart = 10\n",
+    "z = np.repeat(ds_fields.depth[0].values, npart)\n",
+    "lons = np.repeat(31, npart)\n",
+    "lats = np.linspace(-32.5, -30.5, npart)\n",
+    "\n",
+    "pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, lat=lats, lon=lons)\n",
+    "output_file = parcels.ParticleFile(\n",
+    "    store=\"summed_advection_wind.zarr\", outputdt=np.timedelta64(6, \"h\")\n",
+    ")\n",
+    "pset.execute(\n",
+    "    [parcels.kernels.AdvectionRK4],\n",
+    "    runtime=np.timedelta64(5, \"D\"),\n",
+    "    dt=np.timedelta64(1, \"h\"),\n",
+    "    output_file=output_file,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot the resulting particle trajectories overlapped for both cases\n",
+    "summed_advection_wind = xr.open_zarr(\"summed_advection_wind.zarr\")\n",
+    "plt.plot(summed_advection_wind.lon.T, summed_advection_wind.lat.T, \"-\")\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "docs",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.14.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/user_guide/index.md

@@ -29,6 +29,7 @@ examples/tutorial_nemo_curvilinear.ipynb
 examples/tutorial_nemo_3D.ipynb
 examples/tutorial_velocityconversion.ipynb
 examples/tutorial_nestedgrids.ipynb
+examples/tutorial_summingfields.ipynb
 ```
 
 <!-- examples/documentation_indexing.ipynb -->