From 41118e95410c1c1f8cf5c8296e169e42977bd1c8 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 20 May 2026 17:17:42 -0400 Subject: [PATCH 1/3] Add fesom tutorial to documentation This tutorial uses the provided fesom2 channel data for demonstrating the basics of getting started with FESOM2 (unstructured grid) in parcels --- docs/user_guide/examples/tutorial_fesom.ipynb | 204 ++++++++++++++++++ docs/user_guide/index.md | 1 + 2 files changed, 205 insertions(+) create mode 100644 docs/user_guide/examples/tutorial_fesom.ipynb diff --git a/docs/user_guide/examples/tutorial_fesom.ipynb b/docs/user_guide/examples/tutorial_fesom.ipynb new file mode 100644 index 000000000..69e384dcc --- /dev/null +++ b/docs/user_guide/examples/tutorial_fesom.ipynb @@ -0,0 +1,204 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 🖥️ FESOM tutorial\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parcels v4 supports unstructured-grid model output via [uxarray](https://uxarray.readthedocs.io/). This tutorial walks through the minimum steps to advect particles in real [FESOM2](https://fesom.de/) output. The recipe is:\n", + "\n", + "1. Open the FESOM grid and data files with `uxarray`.\n", + "2. Rename FESOM-specific dimensions to Parcels' UGRID conventions with `parcels.convert.fesom_to_ugrid`.\n", + "3. Build a `FieldSet` with `parcels.FieldSet.from_ugrid_conventions`.\n", + "4. Run the simulation as on any structured grid.\n", + "\n", + "If you have not done so already, work through the [quickstart tutorial](../../getting_started/tutorial_quickstart.md) first to get familiar with `ParticleSet`, `Kernel`, and `ParticleFile`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "from pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport uxarray as ux\n\nimport parcels\nimport parcels.tutorial\nfrom parcels.convert import fesom_to_ugrid" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Get the FESOM tutorial dataset\n", + "\n", + "We use a small periodic-channel snapshot from a FESOM2 simulation that ships with Parcels' tutorial data registry. As in the [quickstart](../../getting_started/tutorial_quickstart.md), `parcels.tutorial.open_dataset` downloads the files into a local cache on first use; subsequent calls just return the cached copy.\n", + "\n", + "`uxarray` expects file paths rather than an in-memory dataset, so we trigger the downloads and then point `ux.open_mfdataset` at the cached files:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for name in [\n", + " \"FESOM_periodic_channel/fesom_channel\", # grid description\n", + " \"FESOM_periodic_channel/u.fesom_channel\", # zonal velocity (face-registered)\n", + " \"FESOM_periodic_channel/v.fesom_channel\", # meridional velocity (face-registered)\n", + " \"FESOM_periodic_channel/w.fesom_channel\", # vertical velocity (node-registered)\n", + "]:\n", + " parcels.tutorial.open_dataset(name)\n", + "\n", + "from parcels._datasets.remote import _DATA_HOME\n", + "data_dir = Path(_DATA_HOME) / \"data\" / \"FESOM_periodic_channel\"\n", + "\n", + "grid_path = str(data_dir / \"fesom_channel.nc\")\n", + "data_paths = [\n", + " str(data_dir / \"u.fesom_channel.nc\"),\n", + " str(data_dir / \"v.fesom_channel.nc\"),\n", + " str(data_dir / \"w.fesom_channel.nc\"),\n", + "]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```{note}\n", + "The `FESOM_periodic_channel` dataset is a single-time-step snapshot of an idealised channel (~4.4° × ~18° wide, 40 vertical layers). Working with multi-time FESOM output is identical, except you pass a glob or list of data files to `ux.open_mfdataset`.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Open the data with `uxarray`\n", + "\n", + "`ux.open_mfdataset(grid_path, data_paths)` reads the FESOM grid description and joins the velocity files on its grid. FESOM names its velocity variables `u`, `v`, `w` — we rename them to `U`, `V`, `W` so that `parcels.FieldSet.from_ugrid_conventions` recognises them as the velocity components:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = ux.open_mfdataset(grid_path, data_paths).rename_vars({\"u\": \"U\", \"v\": \"V\", \"w\": \"W\"})\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note the FESOM-specific dimension names: `elem` (number of triangular faces), `nod2` (number of nodes), `nz1` (vertical layer centres), and `nz` (layer interfaces). The horizontal velocities `U` and `V` live on face centres along `nz1`; the vertical velocity `W` lives on nodes along `nz`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Convert to UGRID conventions\n", + "\n", + "Parcels works with a small UGRID-compliant dialect: nodes are `n_node`, faces are `n_face`, vertical centres are `zc`, and vertical interfaces are `zf`. The helper `parcels.convert.fesom_to_ugrid` does this rename in one call:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = fesom_to_ugrid(ds)\n", + "print(\"dims:\", dict(ds.sizes))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Build the `FieldSet`\n", + "\n", + "With UGRID-compliant dimensions in place, `parcels.FieldSet.from_ugrid_conventions` builds the `FieldSet`. It detects `U`, `V`, `W`, assigns a `UxGrid` to each field, and picks an appropriate interpolator based on each variable's coordinate location (face- vs. node-registered, centre vs. interface). Use `mesh=\"spherical\"` so that velocities in m/s are correctly converted to deg/s along the lon/lat coordinates:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh=\"spherical\")\n", + "\n", + "for name, field in fieldset.fields.items():\n", + " interp = getattr(field, \"interp_method\", None)\n", + " interp_name = interp.__name__ if interp is not None else \"-\"\n", + " print(f\"{name:>4s} -> {type(field).__name__:<11s} interp={interp_name}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`U` and `V` get face-registered interpolation, `W` gets node-registered linear interpolation. The combined vector fields `UV` and `UVW` are assembled automatically." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "## Release particles and advect\n\nWe seed particles on a grid of four latitudes spanning the channel and ten longitudes, and integrate for two days with RK4. Because this snapshot has only a single time level, `fieldset.time_interval` is `None` and we omit the `time=` argument so that Parcels treats the flow as constant in time:" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "lon_grid, lat_grid = np.meshgrid(\n np.linspace(0.5, 4.0, 10),\n np.linspace(3.0, 15.0, 4),\n)\nlon = lon_grid.ravel()\nlat = lat_grid.ravel()\nz = np.full(lon.size, 50.0) # release at 50 m depth\n\npset = parcels.ParticleSet(\n fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, z=z,\n)\n\noutput_file = parcels.ParticleFile(\"output-fesom.parquet\", outputdt=np.timedelta64(1, \"h\"))\n\npset.execute(\n [parcels.kernels.AdvectionRK4],\n runtime=np.timedelta64(2, \"D\"),\n dt=np.timedelta64(5, \"m\"),\n output_file=output_file,\n verbose_progress=False,\n)" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "## Plot the result\n\nEach particle trajectory is coloured by observation time so we can follow how the particles drift through the FESOM2 velocity field:" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "df = parcels.read_particlefile(\"output-fesom.parquet\")\n\nfig, ax = plt.subplots(figsize=(10, 5))\nax.scatter(lon, lat, facecolors=\"none\", edgecolors=\"k\", s=60, label=\"release\")\nsc = ax.scatter(df[\"lon\"], df[\"lat\"], c=df[\"time\"].dt.total_seconds(), s=8, cmap=\"viridis\")\nfig.colorbar(sc, ax=ax, label=\"time since release [s]\")\nax.set_xlabel(\"Longitude [deg E]\")\nax.set_ylabel(\"Latitude [deg N]\")\nax.legend(loc=\"upper right\")\nplt.show()" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The particles drift through the channel following the FESOM2 velocity field. From here, the rest of Parcels — custom kernels, sampling fields onto particles, writing your own interpolators — works identically to structured grids. See the [interpolation tutorial](./tutorial_interpolation.ipynb) for the available `Ux*` interpolators and how to write a custom one." + ] + } + ], + "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.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index 3303f1bc6..9e22f1236 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -28,6 +28,7 @@ examples/explanation_grids.md examples/tutorial_nemo.ipynb examples/tutorial_croco_3D.ipynb examples/tutorial_mitgcm.ipynb +examples/tutorial_fesom.ipynb examples/tutorial_velocityconversion.ipynb examples/tutorial_nestedgrids.ipynb examples/tutorial_manipulating_field_data.ipynb From 701ce133495f7e11b5ca5432499a52ea256f97bb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 20 May 2026 21:19:41 +0000 Subject: [PATCH 2/3] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- docs/user_guide/examples/tutorial_fesom.ipynb | 66 +++++++++++++++++-- 1 file changed, 60 insertions(+), 6 deletions(-) diff --git a/docs/user_guide/examples/tutorial_fesom.ipynb b/docs/user_guide/examples/tutorial_fesom.ipynb index 69e384dcc..2246ab8ed 100644 --- a/docs/user_guide/examples/tutorial_fesom.ipynb +++ b/docs/user_guide/examples/tutorial_fesom.ipynb @@ -26,7 +26,17 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "from pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport uxarray as ux\n\nimport parcels\nimport parcels.tutorial\nfrom parcels.convert import fesom_to_ugrid" + "source": [ + "from pathlib import Path\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import uxarray as ux\n", + "\n", + "import parcels\n", + "import parcels.tutorial\n", + "from parcels.convert import fesom_to_ugrid" + ] }, { "cell_type": "markdown", @@ -46,7 +56,7 @@ "outputs": [], "source": [ "for name in [\n", - " \"FESOM_periodic_channel/fesom_channel\", # grid description\n", + " \"FESOM_periodic_channel/fesom_channel\", # grid description\n", " \"FESOM_periodic_channel/u.fesom_channel\", # zonal velocity (face-registered)\n", " \"FESOM_periodic_channel/v.fesom_channel\", # meridional velocity (face-registered)\n", " \"FESOM_periodic_channel/w.fesom_channel\", # vertical velocity (node-registered)\n", @@ -54,6 +64,7 @@ " parcels.tutorial.open_dataset(name)\n", "\n", "from parcels._datasets.remote import _DATA_HOME\n", + "\n", "data_dir = Path(_DATA_HOME) / \"data\" / \"FESOM_periodic_channel\"\n", "\n", "grid_path = str(data_dir / \"fesom_channel.nc\")\n", @@ -88,7 +99,9 @@ "metadata": {}, "outputs": [], "source": [ - "ds = ux.open_mfdataset(grid_path, data_paths).rename_vars({\"u\": \"U\", \"v\": \"V\", \"w\": \"W\"})\n", + "ds = ux.open_mfdataset(grid_path, data_paths).rename_vars(\n", + " {\"u\": \"U\", \"v\": \"V\", \"w\": \"W\"}\n", + ")\n", "ds" ] }, @@ -158,7 +171,35 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "lon_grid, lat_grid = np.meshgrid(\n np.linspace(0.5, 4.0, 10),\n np.linspace(3.0, 15.0, 4),\n)\nlon = lon_grid.ravel()\nlat = lat_grid.ravel()\nz = np.full(lon.size, 50.0) # release at 50 m depth\n\npset = parcels.ParticleSet(\n fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, z=z,\n)\n\noutput_file = parcels.ParticleFile(\"output-fesom.parquet\", outputdt=np.timedelta64(1, \"h\"))\n\npset.execute(\n [parcels.kernels.AdvectionRK4],\n runtime=np.timedelta64(2, \"D\"),\n dt=np.timedelta64(5, \"m\"),\n output_file=output_file,\n verbose_progress=False,\n)" + "source": [ + "lon_grid, lat_grid = np.meshgrid(\n", + " np.linspace(0.5, 4.0, 10),\n", + " np.linspace(3.0, 15.0, 4),\n", + ")\n", + "lon = lon_grid.ravel()\n", + "lat = lat_grid.ravel()\n", + "z = np.full(lon.size, 50.0) # release at 50 m depth\n", + "\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " pclass=parcels.Particle,\n", + " lon=lon,\n", + " lat=lat,\n", + " z=z,\n", + ")\n", + "\n", + "output_file = parcels.ParticleFile(\n", + " \"output-fesom.parquet\", outputdt=np.timedelta64(1, \"h\")\n", + ")\n", + "\n", + "pset.execute(\n", + " [parcels.kernels.AdvectionRK4],\n", + " runtime=np.timedelta64(2, \"D\"),\n", + " dt=np.timedelta64(5, \"m\"),\n", + " output_file=output_file,\n", + " verbose_progress=False,\n", + ")" + ] }, { "cell_type": "markdown", @@ -170,7 +211,20 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "df = parcels.read_particlefile(\"output-fesom.parquet\")\n\nfig, ax = plt.subplots(figsize=(10, 5))\nax.scatter(lon, lat, facecolors=\"none\", edgecolors=\"k\", s=60, label=\"release\")\nsc = ax.scatter(df[\"lon\"], df[\"lat\"], c=df[\"time\"].dt.total_seconds(), s=8, cmap=\"viridis\")\nfig.colorbar(sc, ax=ax, label=\"time since release [s]\")\nax.set_xlabel(\"Longitude [deg E]\")\nax.set_ylabel(\"Latitude [deg N]\")\nax.legend(loc=\"upper right\")\nplt.show()" + "source": [ + "df = parcels.read_particlefile(\"output-fesom.parquet\")\n", + "\n", + "fig, ax = plt.subplots(figsize=(10, 5))\n", + "ax.scatter(lon, lat, facecolors=\"none\", edgecolors=\"k\", s=60, label=\"release\")\n", + "sc = ax.scatter(\n", + " df[\"lon\"], df[\"lat\"], c=df[\"time\"].dt.total_seconds(), s=8, cmap=\"viridis\"\n", + ")\n", + "fig.colorbar(sc, ax=ax, label=\"time since release [s]\")\n", + "ax.set_xlabel(\"Longitude [deg E]\")\n", + "ax.set_ylabel(\"Latitude [deg N]\")\n", + "ax.legend(loc=\"upper right\")\n", + "plt.show()" + ] }, { "cell_type": "markdown", @@ -201,4 +255,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} \ No newline at end of file +} From 8820d215e26f5dd69c1a5c8b700779f2162f118c Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Fri, 29 May 2026 12:08:35 -0400 Subject: [PATCH 3/3] Overlay particle trajectories with velocity field vectors+speed color --- docs/user_guide/examples/tutorial_fesom.ipynb | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/docs/user_guide/examples/tutorial_fesom.ipynb b/docs/user_guide/examples/tutorial_fesom.ipynb index 69e384dcc..69455e98c 100644 --- a/docs/user_guide/examples/tutorial_fesom.ipynb +++ b/docs/user_guide/examples/tutorial_fesom.ipynb @@ -26,7 +26,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": "from pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport uxarray as ux\n\nimport parcels\nimport parcels.tutorial\nfrom parcels.convert import fesom_to_ugrid" + "source": "from pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport matplotlib.tri as mtri\nimport numpy as np\nimport uxarray as ux\n\nimport parcels\nimport parcels.tutorial" }, { "cell_type": "markdown", @@ -113,10 +113,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "ds = fesom_to_ugrid(ds)\n", - "print(\"dims:\", dict(ds.sizes))" - ] + "source": "ds = parcels.convert.fesom_to_ugrid(ds)\nprint(\"dims:\", dict(ds.sizes))" }, { "cell_type": "markdown", @@ -163,14 +160,14 @@ { "cell_type": "markdown", "metadata": {}, - "source": "## Plot the result\n\nEach particle trajectory is coloured by observation time so we can follow how the particles drift through the FESOM2 velocity field:" + "source": "## Plot the velocity field and trajectories\n\nWe plot the particle paths on top of the velocity field they advect through: triangle colour shows the speed at the release depth (≈50 m), black arrows show the velocity at face centres (drawn in lon/lat space, length proportional to speed), grey lines trace each particle's path, and the coloured dots mark the positions over time. Drawing the arrows with `angles=\"xy\"` keeps them aligned with the trajectories, so you can see the particles streak along the fast jets and barely move in the quiet bands between them:" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": "df = parcels.read_particlefile(\"output-fesom.parquet\")\n\nfig, ax = plt.subplots(figsize=(10, 5))\nax.scatter(lon, lat, facecolors=\"none\", edgecolors=\"k\", s=60, label=\"release\")\nsc = ax.scatter(df[\"lon\"], df[\"lat\"], c=df[\"time\"].dt.total_seconds(), s=8, cmap=\"viridis\")\nfig.colorbar(sc, ax=ax, label=\"time since release [s]\")\nax.set_xlabel(\"Longitude [deg E]\")\nax.set_ylabel(\"Latitude [deg N]\")\nax.legend(loc=\"upper right\")\nplt.show()" + "source": "df = parcels.read_particlefile(\"output-fesom.parquet\")\n\ntriang = mtri.Triangulation(\n ds.uxgrid.node_lon.values,\n ds.uxgrid.node_lat.values,\n triangles=ds.uxgrid.face_node_connectivity.values,\n)\n\ndepth_idx = int(np.argmin(np.abs(ds.zc.values - 50.0)))\nU_face = np.asarray(ds[\"U\"].isel(zc=depth_idx)).squeeze()\nV_face = np.asarray(ds[\"V\"].isel(zc=depth_idx)).squeeze()\nspeed = np.hypot(U_face, V_face)\n\nfig, ax = plt.subplots(figsize=(11, 5))\n\n# Background: speed at the release depth.\ntpc = ax.tripcolor(triang, facecolors=speed, shading=\"flat\", cmap=\"Blues\")\nfig.colorbar(tpc, ax=ax, label=\"speed [m/s]\", location=\"left\", shrink=0.85)\n\n# Velocity arrows at face centres. Drawing in lon/lat space (angles/scale_units\n# \"xy\") keeps the arrows aligned with the trajectories; length tracks speed.\nstep = max(1, U_face.size // 400)\nxf = ds.uxgrid.face_lon.values\nyf = ds.uxgrid.face_lat.values\nmax_speed = float(np.nanmax(speed))\nq = ax.quiver(\n xf[::step], yf[::step], U_face[::step], V_face[::step],\n angles=\"xy\", scale_units=\"xy\", scale=max_speed / 0.3,\n color=\"k\", width=0.0018, pivot=\"tail\",\n)\nax.quiverkey(q, 0.86, 1.04, max_speed, f\"{max_speed:.2f} m/s\", labelpos=\"E\", coordinates=\"axes\")\n\n# Particle paths (grey lines) and positions coloured by time.\nfor traj in df.sort(\"time\").partition_by(\"particle_id\"):\n ax.plot(traj[\"lon\"], traj[\"lat\"], color=\"0.4\", linewidth=0.6, alpha=0.7, zorder=2)\nax.scatter(lon, lat, facecolors=\"none\", edgecolors=\"k\", s=60, label=\"release\", zorder=3)\nsc = ax.scatter(\n df[\"lon\"], df[\"lat\"], c=df[\"time\"].dt.total_seconds(), s=6, cmap=\"viridis\", zorder=3,\n)\nfig.colorbar(sc, ax=ax, label=\"time since release [s]\", shrink=0.85)\n\nax.set_xlabel(\"Longitude [deg E]\")\nax.set_ylabel(\"Latitude [deg N]\")\nax.set_title(f\"FESOM2 velocity at z ≈ {float(ds.zc.values[depth_idx]):.1f} m with particle trajectories\")\nax.legend(loc=\"upper right\")\nplt.show()" }, { "cell_type": "markdown",