Skip to content

[DOCS] Delft3D-FM instructions #2992

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
erikvansebille merged 11 commits into from
Oct 20, 2026
Merged

[DOCS] Delft3D-FM instructions #2992

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
8710307
Adding Delft3D-FM instructions to tutorial
erikvansebille Oct 19, 2026
37ea319
Adding remote for Delft3D-fm tutorial
erikvansebille Oct 19, 2026
8b59f3e
Increase figure size
erikvansebille Oct 19, 2026
84e4636
Further cleaning up code
erikvansebille Oct 19, 2026
fcea588
Separate plotting cell
erikvansebille Oct 19, 2026
67515cf
Tidy Delft3D-FM tutorial prose and split conversion cell
VeckoTheGecko Oct 20, 2026
1214db3
Use fan triangulation for Delft3D-FM faces
VeckoTheGecko Oct 20, 2026
da7e3b1
Copy edit
VeckoTheGecko Oct 20, 2026
7df584a
Merge branch 'main' into delft3d-fm_instructions
erikvansebille Oct 20, 2026
bf489f2
Further cleanup notebook
erikvansebille Oct 20, 2026
b11c280
Adding strict=True for zip statement
erikvansebille Oct 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
256 changes: 245 additions & 11 deletions docs/user_guide/examples/tutorial_delft3d.ipynb
View file
Open in desktop

@VeckoTheGecko VeckoTheGecko Oct 20, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Delft3D-fm can also use quadrilateral meshes.

Not quite accurate. Delft 3D FM can also support higher order polygons (see chapter 2 intro at https://content.oss.deltares.nl/delft3d/D-Flow_FM_User_Manual.pdf , or 23.1 on topological conventions "Such a grid cell should contain at least three and at most six corner nodes, representing triangles, quadrangles, pentagons and hexagons, respectively").

The code would have holes in the mesh in the case of these higher order polygons.

I've made an update in 1214db3 accordingly, and to also use fan triangulation.

However, the Parcels results will only be independent of this regridding if you use UxConstantFace... interpolation.

Add "or if your output data is all node-registered" (?)

@VeckoTheGecko VeckoTheGecko Oct 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While reviewing I found that the raw Delft3D-FM tutorial file (delft3d-fm-raw.nc) has inconsistent UGRID metadata, which is why uxarray can't open it directly. Perhaps worth adding a note? (we need to triangular anyway, so opening directly with uxarray isn't an option)

The following was drafted with LLM assistance.

import pooch
import uxarray as ux

path = pooch.retrieve(
    "https://github.com/Parcels-code/parcels-data/raw/main/data/Delft3D_data/delft3d-fm-raw.nc",
    known_hash=None,
)
ux.open_dataset(path, path)
# ValueError: cannot rename 'mesh2d_edge_faces' because it is not a variable or dimension in this dataset

The attributes on the Mesh topology variable don't match the file's contents:

Mesh attribute Value Actually in file
edge_face_connectivity mesh2d_edge_faces variable doesn't exist
node_dimension mesh2d_nNodes nMesh_node
face_dimension mesh2d_nFaces nMesh_face
edge_dimension mesh2d_nEdges nMesh_edge
max_face_nodes_dimension mesh2d_nMax_face_nodes nMaxMesh_face_nodes

It looks like the variables and dimensions were renamed on export (the history attribute mentions FEWS), but the mesh metadata wasn't updated. If I patch the attributes, uxarray reads the grid without problems:

import xarray as xr

ds = xr.open_dataset(path)
attrs = ds["Mesh"].attrs
attrs.pop("edge_face_connectivity")
attrs.update(
    node_dimension="nMesh_node",
    face_dimension="nMesh_face",
    edge_dimension="nMesh_edge",
    max_face_nodes_dimension="nMaxMesh_face_nodes",
)
grid = ux.Grid.from_dataset(ds)  # 7271 faces, 7259 nodes, n_max_face_nodes=4

@erikvansebille erikvansebille Oct 20, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I updated the file to be xgrid-compliant in Parcels-code/parcels-data#15

Original file line number Diff line number Diff line change
Expand Up @@ -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)."
]
},
{
Expand All @@ -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"
Expand All @@ -39,13 +40,246 @@
"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 start with opening the raw file"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4",
"metadata": {},
"outputs": [],
"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"
]
},
{
"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",
"]\n",
"\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",
"ax.set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"id": "7",
"metadata": {},
"source": [
"We now show how to create a regridded version that only contains triangular grid cells."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8",
"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",
" 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",
"\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",
"\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\": ds[\"zc\"]},\n",
" attrs=src.attrs,\n",
" )\n",
"\n",
"\n",
"uxds = ux.UxDataset(\n",
" xr.Dataset(\n",
" {\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",
" ),\n",
" uxgrid=uxgrid,\n",
")\n",
"fieldset = parcels.FieldSet.from_ugrid_conventions(uxds, mesh=\"spherical\")\n",
"fieldset.describe()"
]
},
{
"cell_type": "markdown",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10",
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(8, 6))\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": "11",
"metadata": {},
"source": [
"And can do a simple Parcels simulation with the regridded version of the model output."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12",
"metadata": {},
"outputs": [],
"source": [
"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",
" \"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",
" verbose_progress=False,\n",
")"
]
},
{
"cell_type": "markdown",
"id": "13",
"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": "14",
"metadata": {},
"outputs": [],
"source": [
"df = parcels.read_particlefile(outfile._path)\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",
" ax.plot(traj[\"x\"], traj[\"y\"], \"b.-\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "15",
"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": "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."
]
},
{
"cell_type": "markdown",
"id": "4",
"id": "17",
"metadata": {},
"source": [
"```{note}\n",
Expand All @@ -56,7 +290,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "5",
"id": "18",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -76,7 +310,7 @@
},
{
"cell_type": "markdown",
"id": "6",
"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."
Expand All @@ -85,7 +319,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "7",
"id": "20",
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -123,7 +357,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "8",
"id": "21",
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -168,7 +402,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Parcels:docs (3.14.6)",
"display_name": "Parcels:default (3.14.8)",
"language": "python",
"name": "python3"
},
Expand All @@ -182,7 +416,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.6"
"version": "3.14.8"
},
"mystnb": {
"execution_in_temp": true
Expand Down
2 changes: 2 additions & 0 deletions src/parcels/_datasets/remote.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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)),
Expand Down
Loading