From 87103073310c3bd02b3335e7c3b84bfc43f3b750 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 19 Oct 2026 16:27:54 +0200 Subject: [PATCH 01/10] Adding Delft3D-FM instructions to tutorial --- .../examples/tutorial_delft3d.ipynb | 250 +++++++++++++++++- 1 file changed, 239 insertions(+), 11 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 18f9c9d5b..56f4dc680 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -13,10 +13,7 @@ "id": "1", "metadata": {}, "source": [ - "This tutorial shows how to load in [Delft3D](https://www.deltares.nl/en/software-and-data/products/delft3d-4-suite) model output files (in NetCDF format) into Parcels.\n", - "\n", - "## Structured Grids\n", - "Special about Delft3D is that its structured grid contains NaN values for points that don't exist in the model domain. The hashtable approach in Parcels v4 supports these NaN grid cells out of the box, so that you can use the model output directly without any preprocessing." + "This tutorial shows how to load in [Delft3D](https://www.deltares.nl/en/software-and-data/products/delft3d-4-suite) model output files (in NetCDF format) into Parcels. It treats both the [unstructured grid format (Delft3D-fm)](#unstructured-grids-delft3d-fm) and the [structured grid format (Delft3D-flow)](#structured-grids-delft3d-flow)." ] }, { @@ -28,8 +25,12 @@ "source": [ "import cmocean as cmo\n", "import matplotlib.pyplot as plt\n", + "import matplotlib.tri as mtri\n", "import numpy as np\n", "import pandas as pd\n", + "import uxarray as ux\n", + "import xarray as xr\n", + "from matplotlib.collections import PolyCollection\n", "\n", "import parcels\n", "import parcels.tutorial" @@ -39,13 +40,240 @@ "cell_type": "markdown", "id": "3", "metadata": {}, + "source": [ + "## Unstructured Grids (Delft3D-fm)\n", + "\n", + "Delft3D-fm is the flexible mesh version of Delft3D. It uses an unstructured grid, which is supported in Parcels v4. \n", + "\n", + "```{note}\n", + "There is one major caveat to using Delft3D-fm with Parcels: currently, Parcels v4 can only handle triangular meshes, while Delft3D-fm can also use quadrilateral meshes. Therefore, if you want to use Parcels v4 with Delft3D-fm, you need to first modify your model output to be based on a triangular mesh. Any quadrilateral grid cell can be split into two triangular grid cells, for example by drawing a diagonal line between two opposite corners of the quadrilateral. However, the Parcels results will only be independent of this regridding if you use `UxConstantFace...` interpolation.\n", + "```\n", + "\n", + "The following example shows how to load in a Delft3D-fm model output file and create a `FieldSet` from it. \n", + "\n", + "We first show the raw file that is a combination of triangular and quadrilateral grid cells. This is a small region of the North Sea, just offshore of the city of The Hague." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "raw_ds = xr.open_dataset(\"/users/erik/Downloads/delt3d-fm-raw.nc\")\n", + "\n", + "node_x = raw_ds[\"Mesh_node_x\"].values\n", + "node_y = raw_ds[\"Mesh_node_y\"].values\n", + "raw_face_nodes = raw_ds[\"Mesh_face_nodes\"].values\n", + "\n", + "raw_polygons = [\n", + " [(node_x[i], node_y[i]) for i in face[np.isfinite(face)].astype(np.int64) - 1]\n", + " for face in raw_face_nodes\n", + "]\n", + "\n", + "fig, ax = plt.subplots(figsize=(5, 3))\n", + "ax.add_collection(\n", + " PolyCollection(raw_polygons, edgecolor=\"k\", facecolor=\"none\", linewidth=0.3)\n", + ")\n", + "ax.set_aspect(\"equal\")" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "We now show how to create a regridded version that only contains triangular grid cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "# Use integer-safe connectivity for triangulation logic (1-based with 0 as invalid).\n", + "raw_face_nodes_i = np.where(np.isfinite(raw_face_nodes), raw_face_nodes, 0).astype(\n", + " \"int64\"\n", + ")\n", + "\n", + "tri_list = []\n", + "tri_parent_face = []\n", + "\n", + "for i, face in enumerate(raw_face_nodes_i):\n", + " # UGRID connectivity is typically 1-based with 0/-1 as fill values.\n", + " valid = face[face > 0] - 1\n", + " n_valid = valid.size\n", + "\n", + " if n_valid < 3:\n", + " continue\n", + " if n_valid == 3:\n", + " tri_list.append(valid)\n", + " tri_parent_face.append(i)\n", + " continue\n", + "\n", + " if n_valid == 4:\n", + " a, b, c, d = valid\n", + " diag_ac = (node_x[a] - node_x[c]) ** 2 + (node_y[a] - node_y[c]) ** 2\n", + " diag_bd = (node_x[b] - node_x[d]) ** 2 + (node_y[b] - node_y[d]) ** 2\n", + "\n", + " if diag_ac <= diag_bd:\n", + " tri_list.extend([[a, b, c], [a, c, d]])\n", + " else:\n", + " tri_list.extend([[a, b, d], [b, c, d]])\n", + " tri_parent_face.extend([i, i])\n", + " continue\n", + "\n", + " # Fallback for polygons with >4 nodes: simple fan triangulation.\n", + " for j in range(1, n_valid - 1):\n", + " tri_list.append([valid[0], valid[j], valid[j + 1]])\n", + " tri_parent_face.append(i)\n", + "\n", + "tri_faces = np.asarray(tri_list, dtype=\"int64\")\n", + "tri_parent_face = np.asarray(tri_parent_face, dtype=\"int64\")\n", + "\n", + "zf = [0, 1]\n", + "zc = [0.5]\n", + "\n", + "uxgrid = ux.Grid.from_topology(\n", + " node_lon=raw_ds[\"Mesh_node_x\"],\n", + " node_lat=raw_ds[\"Mesh_node_y\"],\n", + " face_node_connectivity=tri_faces,\n", + " fill_value=-1,\n", + ")\n", + "\n", + "# Fast remap from original faces -> triangular faces using NumPy indexing.\n", + "tri_parent_idx = tri_parent_face.astype(np.intp, copy=False)\n", + "\n", + "\n", + "def to_face_3d(ds, var_name, target_name):\n", + " src = ds[var_name].transpose(\"time\", \"nMesh_face\").astype(\"float32\")\n", + " tri_vals = np.take(src.values, tri_parent_idx, axis=1)\n", + " return xr.DataArray(\n", + " tri_vals[:, np.newaxis, :],\n", + " dims=(\"time\", \"zc\", \"n_face\"),\n", + " coords={\"time\": ds[\"time\"], \"zc\": zc},\n", + " name=target_name,\n", + " attrs=src.attrs,\n", + " )\n", + "\n", + "\n", + "uxds = ux.UxDataset(\n", + " xr.Dataset(\n", + " {\n", + " \"U\": to_face_3d(raw_ds, \"velu\", \"U\"),\n", + " \"V\": to_face_3d(raw_ds, \"velv\", \"V\"),\n", + " },\n", + " coords={\"zf\": (\"zf\", zf), \"zc\": (\"zc\", zc), \"time\": raw_ds[\"time\"]},\n", + " ),\n", + " uxgrid=uxgrid,\n", + ")\n", + "fieldset = parcels.FieldSet.from_ugrid_conventions(uxds, mesh=\"spherical\")\n", + "fieldset.describe()" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "We can also visualise this regridded version of the model output, which is now compatible with Parcels v4. As you can see, all the quadrilateral grid cells have been split into two triangular grid cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(5, 3))\n", + "triang = mtri.Triangulation(\n", + " uxgrid[\"node_lon\"].data,\n", + " uxgrid[\"node_lat\"].data,\n", + " triangles=uxgrid[\"face_node_connectivity\"].data,\n", + ")\n", + "\n", + "ax.triplot(triang, color=\"k\", lw=0.3, alpha=0.5, label=\"Flow mesh\", zorder=2)\n", + "ax.set_aspect(\"equal\")" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "And can do a simple Parcels simulation with the regridded version of the model output." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "X, Y = np.meshgrid(np.linspace(4.02, 4.18, 3), np.linspace(52.05, 52.15, 3))\n", + "pset = parcels.ParticleSet(fieldset, x=X, y=Y)\n", + "\n", + "outfile = parcels.ParticleFile(\n", + " \"tutorial_delft3d_particles.parquet\", outputdt=np.timedelta64(15, \"m\"), mode=\"w\"\n", + ")\n", + "\n", + "pset.execute(\n", + " parcels.kernels.AdvectionRK2,\n", + " dt=np.timedelta64(60, \"s\"),\n", + " endtime=fieldset.time_interval.right,\n", + " output_file=outfile,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "And then we can plot the trajectories of the particles on top of the flow mesh." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "df = parcels.read_particlefile(outfile._path)\n", + "fig, ax = plt.subplots(figsize=(5, 3))\n", + "ax.triplot(triang, color=\"k\", lw=0.3, alpha=0.5, label=\"Flow mesh\", zorder=2)\n", + "ax.set_aspect(\"equal\")\n", + "for traj in df.partition_by(\"particle_id\", maintain_order=True):\n", + " ax.plot(traj[\"x\"], traj[\"y\"], \"b.-\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Structured Grids (Delft3D-flow)\n", + "Special about Delft3D is that its structured grid contains NaN values for points that don't exist in the model domain. The hashtable approach in Parcels v4 supports these NaN grid cells out of the box, so that you can use the model output directly without any preprocessing." + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, "source": [ "The example below is of a small domain in the Port of Rotterdam. We first pick the coorindates and the velocity fields from the dataset, and then convert them to the sgrid conventions. Finally, we create a FieldSet from the converted dataset." ] }, { "cell_type": "markdown", - "id": "4", + "id": "15", "metadata": {}, "source": [ "```{note}\n", @@ -56,7 +284,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -76,7 +304,7 @@ }, { "cell_type": "markdown", - "id": "6", + "id": "17", "metadata": {}, "source": [ "Now we define a grid of a few particles to release in the domain, and run a simple advection simulation. The particles are advected by the Delft3D velocity fields, and we can visualize their trajectories." @@ -85,7 +313,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -123,7 +351,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "19", "metadata": {}, "outputs": [], "source": [ @@ -168,7 +396,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Parcels:docs (3.14.6)", + "display_name": "Parcels:default (3.14.8)", "language": "python", "name": "python3" }, @@ -182,7 +410,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.6" + "version": "3.14.8" }, "mystnb": { "execution_in_temp": true From 37ea3192fe8ba8a9c48546a994565bc6180b0bdf Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 19 Oct 2026 17:11:13 +0200 Subject: [PATCH 02/10] Adding remote for Delft3D-fm tutorial --- docs/user_guide/examples/tutorial_delft3d.ipynb | 2 +- src/parcels/_datasets/remote.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 56f4dc680..d9945de8e 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -61,7 +61,7 @@ "metadata": {}, "outputs": [], "source": [ - "raw_ds = xr.open_dataset(\"/users/erik/Downloads/delt3d-fm-raw.nc\")\n", + "raw_ds = parcels.tutorial.open_dataset(\"Delft3D_data/delft3d-fm-raw\")\n", "\n", "node_x = raw_ds[\"Mesh_node_x\"].values\n", "node_y = raw_ds[\"Mesh_node_y\"].values\n", diff --git a/src/parcels/_datasets/remote.py b/src/parcels/_datasets/remote.py index 3e12bb99d..e19454053 100644 --- a/src/parcels/_datasets/remote.py +++ b/src/parcels/_datasets/remote.py @@ -58,6 +58,7 @@ def _get_data_home() -> Path: # ] + [ "data/Delft3D_data/Rotterdam_tiny.nc", + "data/Delft3D_data/delft3d-fm-raw.nc", ] + [ "data/FESOM_periodic_channel/fesom_channel.nc", @@ -234,6 +235,7 @@ class _Purpose(enum.Enum): ("CopernicusMarine_data_for_Argo_tutorial/data", (_V3Dataset(_ODIE,"data/CopernicusMarine_data_for_Argo_tutorial/cmems_mod_glo_phy-*.nc"), _Purpose.TUTORIAL)), ("CopernicusMarine_data_for_Argo_tutorial/data_wind", (_V3Dataset(_ODIE,"data/CopernicusMarine_data_for_Argo_tutorial/cmems_obs-wind_glo_phy*.nc"), _Purpose.TUTORIAL)), ("Delft3D_data/Rotterdam_tiny", (_V3Dataset(_ODIE,"data/Delft3D_data/Rotterdam_tiny.nc"), _Purpose.TUTORIAL)), + ("Delft3D_data/delft3d-fm-raw", (_V3Dataset(_ODIE,"data/Delft3D_data/delft3d-fm-raw.nc"), _Purpose.TUTORIAL)), ("CopernicusMarine_data_for_stuck_particles_tutorial/data", (_V3Dataset(_ODIE,"data/CopernicusMarine_data_for_stuck_particles_tutorial/cmems_mod_glo_phy_my_0.083deg_P1D-m_NL.nc"), _Purpose.TUTORIAL)), # ("DecayingMovingEddy_data/U", (_V3Dataset(_ODIE,"data/DecayingMovingEddy_data/decaying_moving_eddyU.nc"), _Purpose.TUTORIAL)), # ("DecayingMovingEddy_data/V", (_V3Dataset(_ODIE,"data/DecayingMovingEddy_data/decaying_moving_eddyV.nc"), _Purpose.TUTORIAL)), From 8b59f3e11c2d7436241551a961302d1e4cd58ba7 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 19 Oct 2026 17:11:36 +0200 Subject: [PATCH 03/10] Increase figure size --- docs/user_guide/examples/tutorial_delft3d.ipynb | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index d9945de8e..085a0a809 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -72,7 +72,7 @@ " for face in raw_face_nodes\n", "]\n", "\n", - "fig, ax = plt.subplots(figsize=(5, 3))\n", + "fig, ax = plt.subplots(figsize=(8, 6))\n", "ax.add_collection(\n", " PolyCollection(raw_polygons, edgecolor=\"k\", facecolor=\"none\", linewidth=0.3)\n", ")\n", @@ -189,7 +189,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots(figsize=(5, 3))\n", + "fig, ax = plt.subplots(figsize=(8, 6))\n", "triang = mtri.Triangulation(\n", " uxgrid[\"node_lon\"].data,\n", " uxgrid[\"node_lat\"].data,\n", @@ -215,7 +215,7 @@ "metadata": {}, "outputs": [], "source": [ - "X, Y = np.meshgrid(np.linspace(4.02, 4.18, 3), np.linspace(52.05, 52.15, 3))\n", + "X, Y = np.meshgrid(np.linspace(4.02, 4.18, 5), np.linspace(52.05, 52.15, 5))\n", "pset = parcels.ParticleSet(fieldset, x=X, y=Y)\n", "\n", "outfile = parcels.ParticleFile(\n", @@ -227,6 +227,7 @@ " dt=np.timedelta64(60, \"s\"),\n", " endtime=fieldset.time_interval.right,\n", " output_file=outfile,\n", + " verbose_progress=False,\n", ")" ] }, @@ -246,7 +247,7 @@ "outputs": [], "source": [ "df = parcels.read_particlefile(outfile._path)\n", - "fig, ax = plt.subplots(figsize=(5, 3))\n", + "fig, ax = plt.subplots(figsize=(8, 6))\n", "ax.triplot(triang, color=\"k\", lw=0.3, alpha=0.5, label=\"Flow mesh\", zorder=2)\n", "ax.set_aspect(\"equal\")\n", "for traj in df.partition_by(\"particle_id\", maintain_order=True):\n", From 84e46363ef55cbc548a150fa4369e6753b417ece Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 19 Oct 2026 17:29:09 +0200 Subject: [PATCH 04/10] Further cleaning up code --- .../examples/tutorial_delft3d.ipynb | 49 +++++++------------ 1 file changed, 19 insertions(+), 30 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 085a0a809..7e481cc05 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -63,6 +63,9 @@ "source": [ "raw_ds = parcels.tutorial.open_dataset(\"Delft3D_data/delft3d-fm-raw\")\n", "\n", + "# add depth coorinates to the raw dataset -as this is not present in the raw dataset, but is required for Parcels to create a FieldSet\n", + "raw_ds = raw_ds.assign_coords(zf=[0, 1], zc=[0.5])\n", + "\n", "node_x = raw_ds[\"Mesh_node_x\"].values\n", "node_y = raw_ds[\"Mesh_node_y\"].values\n", "raw_face_nodes = raw_ds[\"Mesh_face_nodes\"].values\n", @@ -99,7 +102,7 @@ " \"int64\"\n", ")\n", "\n", - "tri_list = []\n", + "tri_faces = []\n", "tri_parent_face = []\n", "\n", "for i, face in enumerate(raw_face_nodes_i):\n", @@ -110,32 +113,17 @@ " if n_valid < 3:\n", " continue\n", " if n_valid == 3:\n", - " tri_list.append(valid)\n", + " tri_faces.append(valid)\n", " tri_parent_face.append(i)\n", - " continue\n", - "\n", " if n_valid == 4:\n", " a, b, c, d = valid\n", " diag_ac = (node_x[a] - node_x[c]) ** 2 + (node_y[a] - node_y[c]) ** 2\n", " diag_bd = (node_x[b] - node_x[d]) ** 2 + (node_y[b] - node_y[d]) ** 2\n", - "\n", " if diag_ac <= diag_bd:\n", - " tri_list.extend([[a, b, c], [a, c, d]])\n", + " tri_faces.extend([[a, b, c], [a, c, d]])\n", " else:\n", - " tri_list.extend([[a, b, d], [b, c, d]])\n", + " tri_faces.extend([[a, b, d], [b, c, d]])\n", " tri_parent_face.extend([i, i])\n", - " continue\n", - "\n", - " # Fallback for polygons with >4 nodes: simple fan triangulation.\n", - " for j in range(1, n_valid - 1):\n", - " tri_list.append([valid[0], valid[j], valid[j + 1]])\n", - " tri_parent_face.append(i)\n", - "\n", - "tri_faces = np.asarray(tri_list, dtype=\"int64\")\n", - "tri_parent_face = np.asarray(tri_parent_face, dtype=\"int64\")\n", - "\n", - "zf = [0, 1]\n", - "zc = [0.5]\n", "\n", "uxgrid = ux.Grid.from_topology(\n", " node_lon=raw_ds[\"Mesh_node_x\"],\n", @@ -144,18 +132,15 @@ " fill_value=-1,\n", ")\n", "\n", - "# Fast remap from original faces -> triangular faces using NumPy indexing.\n", - "tri_parent_idx = tri_parent_face.astype(np.intp, copy=False)\n", "\n", - "\n", - "def to_face_3d(ds, var_name, target_name):\n", - " src = ds[var_name].transpose(\"time\", \"nMesh_face\").astype(\"float32\")\n", - " tri_vals = np.take(src.values, tri_parent_idx, axis=1)\n", + "def to_face(ds, var_name):\n", + " \"\"\"Utility function to convert a variable to the new grid\"\"\"\n", + " src = ds[var_name].transpose(\"time\", \"nMesh_face\")\n", + " tri_vals = np.take(src.values, tri_parent_face, axis=1)\n", " return xr.DataArray(\n", " tri_vals[:, np.newaxis, :],\n", " dims=(\"time\", \"zc\", \"n_face\"),\n", - " coords={\"time\": ds[\"time\"], \"zc\": zc},\n", - " name=target_name,\n", + " coords={\"time\": ds[\"time\"], \"zc\": ds[\"zc\"]},\n", " attrs=src.attrs,\n", " )\n", "\n", @@ -163,10 +148,14 @@ "uxds = ux.UxDataset(\n", " xr.Dataset(\n", " {\n", - " \"U\": to_face_3d(raw_ds, \"velu\", \"U\"),\n", - " \"V\": to_face_3d(raw_ds, \"velv\", \"V\"),\n", + " \"U\": to_face(raw_ds, \"velu\"),\n", + " \"V\": to_face(raw_ds, \"velv\"),\n", + " },\n", + " coords={\n", + " \"zf\": (\"zf\", raw_ds[\"zf\"].data),\n", + " \"zc\": (\"zc\", raw_ds[\"zc\"].data),\n", + " \"time\": raw_ds[\"time\"],\n", " },\n", - " coords={\"zf\": (\"zf\", zf), \"zc\": (\"zc\", zc), \"time\": raw_ds[\"time\"]},\n", " ),\n", " uxgrid=uxgrid,\n", ")\n", From fcea588ff7f8dd7f3f1d7d7618a9f9a0a1b005a7 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 19 Oct 2026 17:33:05 +0200 Subject: [PATCH 05/10] Separate plotting cell --- .../examples/tutorial_delft3d.ipynb | 52 ++++++++++++------- 1 file changed, 34 insertions(+), 18 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 7e481cc05..94d05b15e 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -51,7 +51,7 @@ "\n", "The following example shows how to load in a Delft3D-fm model output file and create a `FieldSet` from it. \n", "\n", - "We first show the raw file that is a combination of triangular and quadrilateral grid cells. This is a small region of the North Sea, just offshore of the city of The Hague." + "We start with opening the raw file" ] }, { @@ -68,8 +68,24 @@ "\n", "node_x = raw_ds[\"Mesh_node_x\"].values\n", "node_y = raw_ds[\"Mesh_node_y\"].values\n", - "raw_face_nodes = raw_ds[\"Mesh_face_nodes\"].values\n", - "\n", + "raw_face_nodes = raw_ds[\"Mesh_face_nodes\"].values" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "The figure below shows that the raw file is a combination of triangular and quadrilateral grid cells. This is a small region of the North Sea, just offshore of the city of The Hague." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ "raw_polygons = [\n", " [(node_x[i], node_y[i]) for i in face[np.isfinite(face)].astype(np.int64) - 1]\n", " for face in raw_face_nodes\n", @@ -84,7 +100,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "We now show how to create a regridded version that only contains triangular grid cells." @@ -93,7 +109,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -165,7 +181,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "We can also visualise this regridded version of the model output, which is now compatible with Parcels v4. As you can see, all the quadrilateral grid cells have been split into two triangular grid cells." @@ -174,7 +190,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -191,7 +207,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "And can do a simple Parcels simulation with the regridded version of the model output." @@ -200,7 +216,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -222,7 +238,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "And then we can plot the trajectories of the particles on top of the flow mesh." @@ -231,7 +247,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -246,7 +262,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "## Structured Grids (Delft3D-flow)\n", @@ -255,7 +271,7 @@ }, { "cell_type": "markdown", - "id": "14", + "id": "16", "metadata": {}, "source": [ "The example below is of a small domain in the Port of Rotterdam. We first pick the coorindates and the velocity fields from the dataset, and then convert them to the sgrid conventions. Finally, we create a FieldSet from the converted dataset." @@ -263,7 +279,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "```{note}\n", @@ -274,7 +290,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -294,7 +310,7 @@ }, { "cell_type": "markdown", - "id": "17", + "id": "19", "metadata": {}, "source": [ "Now we define a grid of a few particles to release in the domain, and run a simple advection simulation. The particles are advected by the Delft3D velocity fields, and we can visualize their trajectories." @@ -303,7 +319,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -341,7 +357,7 @@ { "cell_type": "code", "execution_count": null, - "id": "19", + "id": "21", "metadata": {}, "outputs": [], "source": [ From 67515cf9bc4950336d90acf05f28c81dda025835 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 20 Oct 2026 11:26:22 +0200 Subject: [PATCH 06/10] Tidy Delft3D-FM tutorial prose and split conversion cell - Read particle file by filename instead of the private outfile._path - Split the triangulation/UxDataset/FieldSet cell into three steps - Fix typos and wording in the Delft3D tutorial --- .../examples/tutorial_delft3d.ipynb | 85 +++++++++++++------ 1 file changed, 59 insertions(+), 26 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 94d05b15e..93a1b2c4c 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -43,15 +43,15 @@ "source": [ "## Unstructured Grids (Delft3D-fm)\n", "\n", - "Delft3D-fm is the flexible mesh version of Delft3D. It uses an unstructured grid, which is supported in Parcels v4. \n", + "Delft3D-fm is the flexible mesh version of Delft3D. It uses an unstructured grid, which is supported in Parcels v4.\n", "\n", "```{note}\n", "There is one major caveat to using Delft3D-fm with Parcels: currently, Parcels v4 can only handle triangular meshes, while Delft3D-fm can also use quadrilateral meshes. Therefore, if you want to use Parcels v4 with Delft3D-fm, you need to first modify your model output to be based on a triangular mesh. Any quadrilateral grid cell can be split into two triangular grid cells, for example by drawing a diagonal line between two opposite corners of the quadrilateral. However, the Parcels results will only be independent of this regridding if you use `UxConstantFace...` interpolation.\n", "```\n", "\n", - "The following example shows how to load in a Delft3D-fm model output file and create a `FieldSet` from it. \n", + "The following example shows how to load in a Delft3D-fm model output file and create a `FieldSet` from it.\n", "\n", - "We start with opening the raw file" + "We start by opening the raw file." ] }, { @@ -63,7 +63,7 @@ "source": [ "raw_ds = parcels.tutorial.open_dataset(\"Delft3D_data/delft3d-fm-raw\")\n", "\n", - "# add depth coorinates to the raw dataset -as this is not present in the raw dataset, but is required for Parcels to create a FieldSet\n", + "# Add depth coordinates to the raw dataset, as these are not present in the raw dataset but are required for Parcels to create a FieldSet\n", "raw_ds = raw_ds.assign_coords(zf=[0, 1], zc=[0.5])\n", "\n", "node_x = raw_ds[\"Mesh_node_x\"].values\n", @@ -103,7 +103,7 @@ "id": "7", "metadata": {}, "source": [ - "We now show how to create a regridded version that only contains triangular grid cells." + "We now show how to create a regridded version that only contains triangular grid cells. First, we split each quadrilateral grid cell into two triangles, keeping track of which original grid cell each triangle came from." ] }, { @@ -139,8 +139,24 @@ " tri_faces.extend([[a, b, c], [a, c, d]])\n", " else:\n", " tri_faces.extend([[a, b, d], [b, c, d]])\n", - " tri_parent_face.extend([i, i])\n", - "\n", + " tri_parent_face.extend([i, i])" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "Next, we build a `uxarray` grid from these triangles, and copy the velocities of each original grid cell onto its triangles." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ "uxgrid = ux.Grid.from_topology(\n", " node_lon=raw_ds[\"Mesh_node_x\"],\n", " node_lat=raw_ds[\"Mesh_node_y\"],\n", @@ -174,14 +190,31 @@ " },\n", " ),\n", " uxgrid=uxgrid,\n", - ")\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "Finally, we create a `FieldSet` from the regridded dataset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ "fieldset = parcels.FieldSet.from_ugrid_conventions(uxds, mesh=\"spherical\")\n", "fieldset.describe()" ] }, { "cell_type": "markdown", - "id": "9", + "id": "13", "metadata": {}, "source": [ "We can also visualise this regridded version of the model output, which is now compatible with Parcels v4. As you can see, all the quadrilateral grid cells have been split into two triangular grid cells." @@ -190,7 +223,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -207,16 +240,16 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "15", "metadata": {}, "source": [ - "And can do a simple Parcels simulation with the regridded version of the model output." + "We can then run a simple Parcels simulation with the regridded version of the model output." ] }, { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -238,20 +271,20 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "17", "metadata": {}, "source": [ - "And then we can plot the trajectories of the particles on top of the flow mesh." + "Finally, we plot the trajectories of the particles on top of the flow mesh." ] }, { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "18", "metadata": {}, "outputs": [], "source": [ - "df = parcels.read_particlefile(outfile._path)\n", + "df = parcels.read_particlefile(\"tutorial_delft3d_particles.parquet\")\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "ax.triplot(triang, color=\"k\", lw=0.3, alpha=0.5, label=\"Flow mesh\", zorder=2)\n", "ax.set_aspect(\"equal\")\n", @@ -262,7 +295,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "19", "metadata": {}, "source": [ "## Structured Grids (Delft3D-flow)\n", @@ -271,15 +304,15 @@ }, { "cell_type": "markdown", - "id": "16", + "id": "20", "metadata": {}, "source": [ - "The example below is of a small domain in the Port of Rotterdam. We first pick the coorindates and the velocity fields from the dataset, and then convert them to the sgrid conventions. Finally, we create a FieldSet from the converted dataset." + "The example below is of a small domain in the Port of Rotterdam. We first pick the coordinates and the velocity fields from the dataset, and then convert them to the sgrid conventions. Finally, we create a FieldSet from the converted dataset." ] }, { "cell_type": "markdown", - "id": "17", + "id": "21", "metadata": {}, "source": [ "```{note}\n", @@ -290,7 +323,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -310,7 +343,7 @@ }, { "cell_type": "markdown", - "id": "19", + "id": "23", "metadata": {}, "source": [ "Now we define a grid of a few particles to release in the domain, and run a simple advection simulation. The particles are advected by the Delft3D velocity fields, and we can visualize their trajectories." @@ -319,7 +352,7 @@ { "cell_type": "code", "execution_count": null, - "id": "20", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -330,7 +363,7 @@ "z = np.full(points_x.shape, 0.06) # Depth of 0.06 m\n", "pset = parcels.ParticleSet(fieldset, x=points_x, y=points_y, z=z)\n", "\n", - "# St up an output file to save particle trajectories\n", + "# Set up an output file to save particle trajectories\n", "outputfile = parcels.ParticleFile(\n", " \"Delft3D_structured.parquet\",\n", " outputdt=np.timedelta64(60, \"s\"),\n", @@ -357,7 +390,7 @@ { "cell_type": "code", "execution_count": null, - "id": "21", + "id": "25", "metadata": {}, "outputs": [], "source": [ From 1214db3622da2b02e71260dc0094f1f363e09c60 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 20 Oct 2026 11:35:26 +0200 Subject: [PATCH 07/10] Use fan triangulation for Delft3D-FM faces The previous loop only handled triangles and quads, silently dropping pentagons and hexagons (which D-Flow FM supports). Also parse the face connectivity once, using the start_index attribute instead of a hard-coded 1-based offset. --- .../examples/tutorial_delft3d.ipynb | 46 +++++++------------ 1 file changed, 16 insertions(+), 30 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 93a1b2c4c..8959794e2 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -46,7 +46,7 @@ "Delft3D-fm is the flexible mesh version of Delft3D. It uses an unstructured grid, which is supported in Parcels v4.\n", "\n", "```{note}\n", - "There is one major caveat to using Delft3D-fm with Parcels: currently, Parcels v4 can only handle triangular meshes, while Delft3D-fm can also use quadrilateral meshes. Therefore, if you want to use Parcels v4 with Delft3D-fm, you need to first modify your model output to be based on a triangular mesh. Any quadrilateral grid cell can be split into two triangular grid cells, for example by drawing a diagonal line between two opposite corners of the quadrilateral. However, the Parcels results will only be independent of this regridding if you use `UxConstantFace...` interpolation.\n", + "There is one major caveat to using Delft3D-fm with Parcels: currently, Parcels v4 can only handle triangular meshes, while Delft3D-fm can also use quadrilateral, pentagonal and hexagonal grid cells. Therefore, if you want to use Parcels v4 with Delft3D-fm, you need to first modify your model output to be based on a triangular mesh. Any convex grid cell with $n$ corners can be split into $n-2$ triangular grid cells, for example by drawing lines from one corner to all other corners (a so-called fan triangulation). However, the Parcels results will only be independent of this regridding if you use `UxConstantFace...` interpolation.\n", "```\n", "\n", "The following example shows how to load in a Delft3D-fm model output file and create a `FieldSet` from it.\n", @@ -68,7 +68,12 @@ "\n", "node_x = raw_ds[\"Mesh_node_x\"].values\n", "node_y = raw_ds[\"Mesh_node_y\"].values\n", - "raw_face_nodes = raw_ds[\"Mesh_face_nodes\"].values" + "# Get the node indices of each face, removing the fill values (NaN) and converting to 0-based indexing\n", + "start_index = int(raw_ds[\"Mesh_face_nodes\"].attrs.get(\"start_index\", 0))\n", + "faces = [\n", + " face[np.isfinite(face)].astype(np.int64) - start_index\n", + " for face in raw_ds[\"Mesh_face_nodes\"].values\n", + "]" ] }, { @@ -86,10 +91,7 @@ "metadata": {}, "outputs": [], "source": [ - "raw_polygons = [\n", - " [(node_x[i], node_y[i]) for i in face[np.isfinite(face)].astype(np.int64) - 1]\n", - " for face in raw_face_nodes\n", - "]\n", + "raw_polygons = [list(zip(node_x[face], node_y[face])) for face in faces]\n", "\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "ax.add_collection(\n", @@ -103,7 +105,7 @@ "id": "7", "metadata": {}, "source": [ - "We now show how to create a regridded version that only contains triangular grid cells. First, we split each quadrilateral grid cell into two triangles, keeping track of which original grid cell each triangle came from." + "We now show how to create a regridded version that only contains triangular grid cells. First, we split each grid cell into triangles using a fan triangulation, keeping track of which original grid cell each triangle came from." ] }, { @@ -113,33 +115,17 @@ "metadata": {}, "outputs": [], "source": [ - "# Use integer-safe connectivity for triangulation logic (1-based with 0 as invalid).\n", - "raw_face_nodes_i = np.where(np.isfinite(raw_face_nodes), raw_face_nodes, 0).astype(\n", - " \"int64\"\n", - ")\n", - "\n", "tri_faces = []\n", "tri_parent_face = []\n", "\n", - "for i, face in enumerate(raw_face_nodes_i):\n", - " # UGRID connectivity is typically 1-based with 0/-1 as fill values.\n", - " valid = face[face > 0] - 1\n", - " n_valid = valid.size\n", - "\n", - " if n_valid < 3:\n", - " continue\n", - " if n_valid == 3:\n", - " tri_faces.append(valid)\n", + "for i, face in enumerate(faces):\n", + " # Fan triangulation: connect the first node to each pair of consecutive other nodes.\n", + " # This is valid for convex faces, and keeps the node ordering (counterclockwise) of the original face\n", + " for k in range(1, len(face) - 1):\n", + " tri_faces.append([face[0], face[k], face[k + 1]])\n", " tri_parent_face.append(i)\n", - " if n_valid == 4:\n", - " a, b, c, d = valid\n", - " diag_ac = (node_x[a] - node_x[c]) ** 2 + (node_y[a] - node_y[c]) ** 2\n", - " diag_bd = (node_x[b] - node_x[d]) ** 2 + (node_y[b] - node_y[d]) ** 2\n", - " if diag_ac <= diag_bd:\n", - " tri_faces.extend([[a, b, c], [a, c, d]])\n", - " else:\n", - " tri_faces.extend([[a, b, d], [b, c, d]])\n", - " tri_parent_face.extend([i, i])" + "\n", + "tri_faces = np.asarray(tri_faces)" ] }, { From da7e3b156d4db38accc17fe3d8473443a4cbd5f9 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 20 Oct 2026 11:43:19 +0200 Subject: [PATCH 08/10] Copy edit --- docs/user_guide/examples/tutorial_delft3d.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 8959794e2..7b07995b8 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -46,7 +46,7 @@ "Delft3D-fm is the flexible mesh version of Delft3D. It uses an unstructured grid, which is supported in Parcels v4.\n", "\n", "```{note}\n", - "There is one major caveat to using Delft3D-fm with Parcels: currently, Parcels v4 can only handle triangular meshes, while Delft3D-fm can also use quadrilateral, pentagonal and hexagonal grid cells. Therefore, if you want to use Parcels v4 with Delft3D-fm, you need to first modify your model output to be based on a triangular mesh. Any convex grid cell with $n$ corners can be split into $n-2$ triangular grid cells, for example by drawing lines from one corner to all other corners (a so-called fan triangulation). However, the Parcels results will only be independent of this regridding if you use `UxConstantFace...` interpolation.\n", + "There is one major caveat to using Delft3D-fm with Parcels: currently, Parcels v4 can only handle triangular meshes, while Delft3D-fm can also use quadrilateral, pentagonal and hexagonal grid cells. Therefore, if you want to use Parcels v4 with Delft3D-fm, you need to first modify your model output to be based on a triangular mesh. Any convex grid cell with $n$ corners can be split into $n-2$ triangular grid cells, for example by drawing lines from one corner to all other corners (i.e., a [fan triangulation](https://en.wikipedia.org/wiki/Fan_triangulation)). However, the Parcels results will only be independent of this regridding if you use `UxConstantFace...` interpolation.\n", "```\n", "\n", "The following example shows how to load in a Delft3D-fm model output file and create a `FieldSet` from it.\n", @@ -421,7 +421,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Parcels:default (3.14.8)", + "display_name": "test-notebooks", "language": "python", "name": "python3" }, @@ -435,7 +435,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.8" + "version": "3.14.7" }, "mystnb": { "execution_in_temp": true From bf489f2cab718f3d369f38fdeb88ad5e2aee54fb Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 20 Oct 2026 16:43:18 +0200 Subject: [PATCH 09/10] Further cleanup notebook --- docs/user_guide/examples/tutorial_delft3d.ipynb | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 7b07995b8..a2de6e355 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -123,9 +123,7 @@ " # This is valid for convex faces, and keeps the node ordering (counterclockwise) of the original face\n", " for k in range(1, len(face) - 1):\n", " tri_faces.append([face[0], face[k], face[k + 1]])\n", - " tri_parent_face.append(i)\n", - "\n", - "tri_faces = np.asarray(tri_faces)" + " tri_parent_face.append(i)" ] }, { @@ -421,7 +419,7 @@ ], "metadata": { "kernelspec": { - "display_name": "test-notebooks", + "display_name": "Parcels:default (3.14.8)", "language": "python", "name": "python3" }, @@ -435,7 +433,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.7" + "version": "3.14.8" }, "mystnb": { "execution_in_temp": true From b11c280151b17583f33fdad364eee9f71f76a59a Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 20 Oct 2026 16:54:08 +0200 Subject: [PATCH 10/10] Adding strict=True for zip statement --- docs/user_guide/examples/tutorial_delft3d.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index a2de6e355..a83329f92 100644 --- a/docs/user_guide/examples/tutorial_delft3d.ipynb +++ b/docs/user_guide/examples/tutorial_delft3d.ipynb @@ -91,7 +91,7 @@ "metadata": {}, "outputs": [], "source": [ - "raw_polygons = [list(zip(node_x[face], node_y[face])) for face in faces]\n", + "raw_polygons = [list(zip(node_x[face], node_y[face], strict=True)) for face in faces]\n", "\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "ax.add_collection(\n",