🖥️ Delft3D tutorial

This tutorial shows how to load in Delft3D model output files (in NetCDF format) into Parcels. It treats both the unstructured grid format (Delft3D-fm) and the structured grid format (Delft3D-flow).

import cmocean as cmo
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np
import pandas as pd
import uxarray as ux
import xarray as xr
from matplotlib.collections import PolyCollection

import parcels
import parcels.tutorial

Unstructured Grids (Delft3D-fm)

Delft3D-fm is the flexible mesh version of Delft3D. It uses an unstructured grid, which is supported in Parcels v4.

Note

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). However, the Parcels results will only be independent of this regridding if you use UxConstantFace... interpolation.

The following example shows how to load in a Delft3D-fm model output file and create a FieldSet from it.

We start by opening the raw file.

raw_ds = parcels.tutorial.open_dataset("Delft3D_data/delft3d-fm-raw")

# 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
raw_ds = raw_ds.assign_coords(zf=[0, 1], zc=[0.5])

node_x = raw_ds["Mesh_node_x"].values
node_y = raw_ds["Mesh_node_y"].values
# Get the node indices of each face, removing the fill values (NaN) and converting to 0-based indexing
start_index = int(raw_ds["Mesh_face_nodes"].attrs.get("start_index", 0))
faces = [
    face[np.isfinite(face)].astype(np.int64) - start_index
    for face in raw_ds["Mesh_face_nodes"].values
]

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.

raw_polygons = [list(zip(node_x[face], node_y[face], strict=True)) for face in faces]

fig, ax = plt.subplots(figsize=(8, 6))
ax.add_collection(
    PolyCollection(raw_polygons, edgecolor="k", facecolor="none", linewidth=0.3)
)
ax.set_aspect("equal")

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.

tri_faces = []
tri_parent_face = []

for i, face in enumerate(faces):
    # Fan triangulation: connect the first node to each pair of consecutive other nodes.
    # This is valid for convex faces, and keeps the node ordering (counterclockwise) of the original face
    for k in range(1, len(face) - 1):
        tri_faces.append([face[0], face[k], face[k + 1]])
        tri_parent_face.append(i)

Next, we build a uxarray grid from these triangles, and copy the velocities of each original grid cell onto its triangles.

uxgrid = ux.Grid.from_topology(
    node_lon=raw_ds["Mesh_node_x"],
    node_lat=raw_ds["Mesh_node_y"],
    face_node_connectivity=tri_faces,
    fill_value=-1,
)


def to_face(ds, var_name):
    """Utility function to convert a variable to the new grid"""
    src = ds[var_name].transpose("time", "nMesh_face")
    tri_vals = np.take(src.values, tri_parent_face, axis=1)
    return xr.DataArray(
        tri_vals[:, np.newaxis, :],
        dims=("time", "zc", "n_face"),
        coords={"time": ds["time"], "zc": ds["zc"]},
        attrs=src.attrs,
    )


uxds = ux.UxDataset(
    xr.Dataset(
        {
            "U": to_face(raw_ds, "velu"),
            "V": to_face(raw_ds, "velv"),
        },
        coords={
            "zf": ("zf", raw_ds["zf"].data),
            "zc": ("zc", raw_ds["zc"].data),
            "time": raw_ds["time"],
        },
    ),
    uxgrid=uxgrid,
)

Finally, we create a FieldSet from the regridded dataset.

fieldset = parcels.FieldSet.from_ugrid_conventions(uxds)
fieldset.describe()

| Name   | Type        |   Grid number | Interp method / value         | Parcels backend   |
|:-------|:------------|--------------:|:------------------------------|:------------------|
| U      | Field       |             0 | UxConstantFaceConstantZC(...) | NumPy             |
| V      | Field       |             0 | UxConstantFaceConstantZC(...) | NumPy             |
| UV     | VectorField |             0 | Ux_Velocity(...)              | -                 |

mesh: SphericalMesh(radius=6366707.019493707)
time interval: (np.datetime64('2025-10-31T21:00:00.000000000'), np.datetime64('2025-10-31T23:00:00.000000000'))

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.

fig, ax = plt.subplots(figsize=(8, 6))
triang = mtri.Triangulation(
    uxgrid["node_lon"].data,
    uxgrid["node_lat"].data,
    triangles=uxgrid["face_node_connectivity"].data,
)

ax.triplot(triang, color="k", lw=0.3, alpha=0.5, label="Flow mesh", zorder=2)
ax.set_aspect("equal")

We can then run a simple Parcels simulation with the regridded version of the model output.

X, Y = np.meshgrid(np.linspace(4.02, 4.18, 5), np.linspace(52.05, 52.15, 5))
pset = parcels.ParticleSet(fieldset, x=X, y=Y)

outfile = parcels.ParticleFile(
    "tutorial_delft3d_particles.parquet", outputdt=np.timedelta64(15, "m"), mode="w"
)

pset.execute(
    parcels.kernels.AdvectionRK2,
    dt=np.timedelta64(60, "s"),
    endtime=fieldset.time_interval.right,
    output_file=outfile,
    verbose_progress=False,
)

INFO: Output files are stored in tutorial_delft3d_particles.parquet

Finally, we plot the trajectories of the particles on top of the flow mesh.

df = parcels.read_particlefile("tutorial_delft3d_particles.parquet")
fig, ax = plt.subplots(figsize=(8, 6))
ax.triplot(triang, color="k", lw=0.3, alpha=0.5, label="Flow mesh", zorder=2)
ax.set_aspect("equal")
for traj in df.partition_by("particle_id", maintain_order=True):
    ax.plot(traj["x"], traj["y"], "b.-")
plt.show()

Structured Grids (Delft3D-flow)

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.

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.

Note

While Delft3D provides velocities on a CGrid, we currently can’t use the parcels.interpolators.CGrid_Velocity interpolator on the velocity fields. This is because the velocities have been rotated to east and north in the output - but the Parcels interpolator expects them in the M and N (along-grid) directions. Therefore, we use the parcels.interpolators.XFreeslip interpolator instead.

ds = parcels.tutorial.open_dataset("Delft3D_data/Rotterdam_tiny")
coords = ds[["XZETA", "YZETA", "SIGMA_C"]]
ds_fset = parcels.convert.delft3d_to_sgrid(
    fields={"U": ds["VELU"], "V": ds["VELV"]}, coords=coords
)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)

# Set the interpolation method for the UV field to XFreeslip (see note above)
fieldset.UV.interp_method = parcels.interpolators.XFreeslip()

fieldset = fieldset.to_windowed_arrays()
fieldset.describe()

| Name   | Type        |   Grid number | Interp method / value   | Parcels backend   |
|:-------|:------------|--------------:|:------------------------|:------------------|
| U      | Field       |             0 | XLinear(...)            | WindowedArray     |
| V      | Field       |             0 | XLinear(...)            | WindowedArray     |
| UV     | VectorField |             0 | XFreeslip(...)          | -                 |

mesh: FlatMesh()
time interval: (np.datetime64('2022-06-21T16:30:00.000000000'), np.datetime64('2022-06-21T21:00:00.000000000'))

/tmp/ipykernel_2657/2007992851.py:3: UserWarning: The delft3d_to_sgrid function is experimental and may not work for all Delft3D datasets. Furthermore, we are not entirely confident that the SGrid layout for Delft3D is implemented correctly. Please report any issues to the Parcels GitHub repository.
  ds_fset = parcels.convert.delft3d_to_sgrid(

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.

# Define a set of starting points
points_x, points_y = np.meshgrid(
    np.linspace(93000, 93300, 5), np.linspace(436300, 436600, 5)
)
z = np.full(points_x.shape, 0.06)  # Depth of 0.06 m
pset = parcels.ParticleSet(fieldset, x=points_x, y=points_y, z=z)

# Set up an output file to save particle trajectories
outputfile = parcels.ParticleFile(
    "Delft3D_structured.parquet",
    outputdt=np.timedelta64(60, "s"),
    mode="w",
)


# Define a custom error handling kernel to delete particles on any error
def DeleteOnAnyError(particles, fieldset):
    any_error = particles.state >= 50  # This captures all Errors
    particles[any_error].state = parcels.StatusCode.Delete


# Run the particle set with the advection kernel and the custom error handling kernel
pset.execute(
    [parcels.kernels.AdvectionRK2, DeleteOnAnyError],
    endtime=fieldset.time_interval.right,
    dt=np.timedelta64(60, "s"),
    output_file=outputfile,
    verbose_progress=False,
)

INFO: Output files are stored in Delft3D_structured.parquet

df = parcels.read_particlefile("Delft3D_structured.parquet")

fig, ax = plt.subplots(figsize=(8, 6))

# Plot background velocity field
time_idx = 0
layer_idx = 0
speed = (
    np.sqrt(ds["VELV"] ** 2 + ds["VELU"] ** 2)
    .isel(TIME=time_idx, LAYER=layer_idx)
    .compute()
)
tpc = plt.pcolor(
    speed.XZETA.values,
    speed.YZETA.values,
    speed,
    edgecolors="black",
    cmap="cmo.speed",
    vmin=0,
    vmax=np.max(np.abs(speed)),
)
time_str = pd.to_datetime(speed.TIME.values).strftime("%Y-%m-%d %H:%M:%S")
fig.colorbar(
    tpc,
    ax=ax,
    label=f"Flow speed on {time_str} and {speed.SIGMA_C.values:.2f}m depth [m/s]",
)

for traj in df.partition_by("particle_id"):
    ax.plot(traj["x"], traj["y"], "b", alpha=0.8)

ax.plot(points_x, points_y, "m.", label="initial particle positions")
ax.set_aspect("equal", adjustable="box")
ax.set_xlabel(ds.XZETA.long_name)
ax.set_ylabel(ds.YZETA.long_name)
plt.show()

/tmp/ipykernel_2657/1231125591.py:13: UserWarning: The input coordinates to pcolor are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolor.
  tpc = plt.pcolor(