Add convert_delft3d_to_sgrid

docs/user_guide/examples/tutorial_delft3d.ipynb

@@ -0,0 +1,190 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# 🖥️ Delft3D tutorial"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1",
+   "metadata": {},
+   "source": [
+    "This tutorial shows how to load in files from the [Delft3D](https://www.deltares.nl/en/software/delft3d-4-suite/) model 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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cmocean as cmo\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "import parcels\n",
+    "import parcels.tutorial"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3",
+   "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",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "While Delft3D provides velocities on a CGrid, we currently can't use the `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 `XFreeslip` interpolator instead.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = parcels.tutorial.open_dataset(\"Delft3D_data/Rotterdam_tiny\")\n",
+    "coords = ds[[\"XZETA\", \"YZETA\", \"SIGMA_C\"]]\n",
+    "ds_fset = parcels.convert.delft3d_to_sgrid(\n",
+    "    fields={\"U\": ds[\"VELU\"], \"V\": ds[\"VELV\"]}, coords=coords\n",
+    ")\n",
+    "fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)\n",
+    "\n",
+    "# Set the interpolation method for the UV field to XFreeslip (see note above)\n",
+    "fieldset.UV.interp_method = parcels.interpolators.XFreeslip()\n",
+    "\n",
+    "fieldset = fieldset.to_windowed_arrays()\n",
+    "fieldset.describe()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Define a set of starting points\n",
+    "points_x, points_y = np.meshgrid(\n",
+    "    np.linspace(93000, 93300, 5), np.linspace(436300, 436600, 5)\n",
+    ")\n",
+    "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",
+    "outputfile = parcels.ParticleFile(\n",
+    "    \"Delft3D_structured.parquet\",\n",
+    "    outputdt=np.timedelta64(60, \"s\"),\n",
+    "    mode=\"w\",\n",
+    ")\n",
+    "\n",
+    "\n",
+    "# Define a custom error handling kernel to delete particles on any error\n",
+    "def DeleteOnAnyError(particles, fieldset):\n",
+    "    any_error = particles.state >= 50  # This captures all Errors\n",
+    "    particles[any_error].state = parcels.StatusCode.Delete\n",
+    "\n",
+    "\n",
+    "# Run the particle set with the advection kernel and the custom error handling kernel\n",
+    "pset.execute(\n",
+    "    [parcels.kernels.AdvectionRK2, DeleteOnAnyError],\n",
+    "    endtime=fieldset.time_interval.right,\n",
+    "    dt=np.timedelta64(60, \"s\"),\n",
+    "    output_file=outputfile,\n",
+    "    verbose_progress=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = parcels.read_particlefile(\"Delft3D_structured.parquet\")\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(8, 6))\n",
+    "\n",
+    "# Plot background velocity field\n",
+    "time_idx = 0\n",
+    "layer_idx = 0\n",
+    "speed = (\n",
+    "    np.sqrt(ds[\"VELV\"] ** 2 + ds[\"VELU\"] ** 2)\n",
+    "    .isel(TIME=time_idx, LAYER=layer_idx)\n",
+    "    .compute()\n",
+    ")\n",
+    "tpc = plt.pcolor(\n",
+    "    speed.XZETA.values,\n",
+    "    speed.YZETA.values,\n",
+    "    speed,\n",
+    "    edgecolors=\"black\",\n",
+    "    cmap=\"cmo.speed\",\n",
+    "    vmin=0,\n",
+    "    vmax=np.max(np.abs(speed)),\n",
+    ")\n",
+    "time_str = pd.to_datetime(speed.TIME.values).strftime(\"%Y-%m-%d %H:%M:%S\")\n",
+    "fig.colorbar(\n",
+    "    tpc,\n",
+    "    ax=ax,\n",
+    "    label=f\"Flow speed on {time_str} and {speed.SIGMA_C.values:.2f}m depth [m/s]\",\n",
+    ")\n",
+    "\n",
+    "for traj in df.partition_by(\"particle_id\"):\n",
+    "    ax.plot(traj[\"x\"], traj[\"y\"], \"b\", alpha=0.8)\n",
+    "\n",
+    "ax.plot(points_x, points_y, \"m.\", label=\"initial particle positions\")\n",
+    "ax.set_aspect(\"equal\", adjustable=\"box\")\n",
+    "ax.set_xlabel(ds.XZETA.long_name)\n",
+    "ax.set_ylabel(ds.YZETA.long_name)\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Parcels:docs (3.14.6)",
+   "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.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/user_guide/index.md

@@ -40,6 +40,7 @@ examples/explanation_grids.md
 examples/tutorial_nemo.ipynb
 examples/tutorial_croco_3D.ipynb
 examples/tutorial_mitgcm.ipynb
+examples/tutorial_delft3d.ipynb
 examples/tutorial_fesom.ipynb
 examples/tutorial_schism.ipynb
 examples/tutorial_velocityconversion.ipynb

pixi.toml

@@ -93,6 +93,7 @@ trajan = "*"
 matplotlib-base = ">=2.0.2"
 gsw = "*"
 py-triangle = "*"
+cmocean = ">=4.0.3,<5"
 
 [feature.devtools.dependencies]
 pdbpp = "*"

src/parcels/_core/xgrid.py

@@ -93,8 +93,9 @@ def _transpose_xfield_data_to_tzyx(da: xr.DataArray, sgrid_metadata: sgrid.SGrid
 
     # All dimensions must be associated with an axis in the grid
     if set(dim_to_axis) != set(da.dims):
+        dims_not_on_grid = set(da.dims) - set(dim_to_axis)
         raise ValueError(
-            f"DataArray {da.name!r} with dims {da.dims} has dimensions that are not associated with a direction on the provided grid."
+            f"DataArray {da.name!r} with dims {da.dims} has dimensions {dims_not_on_grid} that are not associated with a direction on the provided grid."
         )
 
     axes_not_in_field = set(_FIELD_DATA_ORDERING).difference(set(dim_to_axis.values()))

src/parcels/_datasets/remote.py

@@ -55,6 +55,9 @@ def _get_data_home() -> Path:
     #     "data/DecayingMovingEddy_data/decaying_moving_eddyU.nc",
     #     "data/DecayingMovingEddy_data/decaying_moving_eddyV.nc",
     # ]
+    + [
+        "data/Delft3D_data/Rotterdam_tiny.nc",
+    ]
     + [
         "data/FESOM_periodic_channel/fesom_channel.nc",
         "data/FESOM_periodic_channel/u.fesom_channel.nc",
@@ -221,6 +224,7 @@ class _Purpose(enum.Enum):
     # ("Peninsula_data/T", (_V3Dataset(_ODIE,"data/Peninsula_data/peninsulaT.nc"), _Purpose.TUTORIAL)),
     # ("GlobCurrent_example_data/data", (_V3Dataset(_ODIE,"data/GlobCurrent_example_data/*000000-GLOBCURRENT-L4-CUReul_hs-ALT_SUM-v02.0-fv01.0.nc", pre_decode_cf_callable=patch_dataset_v4_compat), _Purpose.TUTORIAL)),
     ("CopernicusMarine_data_for_Argo_tutorial/data", (_V3Dataset(_ODIE,"data/CopernicusMarine_data_for_Argo_tutorial/cmems_mod_glo_phy-*.nc"), _Purpose.TUTORIAL)),
+    ("Delft3D_data/Rotterdam_tiny", (_V3Dataset(_ODIE,"data/Delft3D_data/Rotterdam_tiny.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)),

src/parcels/convert.py

@@ -102,6 +102,23 @@ class _Status(enum.Enum):
     "T": "time",
 }
 
+_DELFT3D_X_EXPECTED_COORDS: list[tuple[str, _Status]] = [(name, _Status.REQUIRED) for name in ["XZETA", "YZETA"]]
+
+_DELFT3D_X_VARNAMES_MAPPING: dict[str, str] = {
+    "XZETA": "lon",
+    "YZETA": "lat",
+    "SIGMA_C": "depth",
+    "TIME": "time",
+}
+
+_DELFT3D_X_AXIS_VARNAMES: dict[str, XgcmAxisDirection] = {
+    "M": "X",
+    "N": "Y",
+    "LAYER": "Z",
+    "time": "T",
+}
+
+
 _CROCO_EXPECTED_COORDS: list[tuple[str, _Status]] = [
     (name, _Status.REQUIRED) for name in ["x_rho", "y_rho", "s_w", "time"]
 ]
@@ -576,6 +593,65 @@ def copernicusmarine_to_sgrid(
     return ds
 
 
+def delft3d_to_sgrid(*, fields: dict[str, xr.Dataset | xr.DataArray], coords: xr.Dataset) -> xr.Dataset:
+    """Create an sgrid-compliant xarray.Dataset from a dataset of structured-grid Delft3D netcdf files.
+
+    Parameters
+    ----------
+    fields : dict[str, xr.Dataset | xr.DataArray]
+        Dictionary of xarray.DataArray objects as obtained from a set of structured-grid Delft3D netcdf files.
+    coords : xarray.Dataset
+        xarray.Dataset containing coordinate variables.
+
+    Returns
+    -------
+    xarray.Dataset
+        Dataset object following SGRID conventions to be (optionally) modified and passed to a FieldSet constructor.
+
+    """
+    warnings.warn(
+        "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.",
+        UserWarning,
+        stacklevel=2,
+    )
+
+    fields = fields.copy()
+    coords = _pick_expected_coords(coords, _DELFT3D_X_EXPECTED_COORDS)
+
+    for name, field_da in fields.items():
+        if isinstance(field_da, xr.Dataset):
+            field_da = field_da[name]
+            # TODO: logging message, warn if multiple fields are in this dataset
+        else:
+            field_da = field_da.rename(name)
+        fields[name] = field_da
+
+    ds = xr.merge(list(fields.values()) + [coords], compat="override")
+
+    ds = _maybe_rename_variables(ds, _DELFT3D_X_VARNAMES_MAPPING)
+    ds = _set_coords(ds, _DELFT3D_X_AXIS_VARNAMES.keys())
+    ds = _set_axis_attrs(ds, _DELFT3D_X_AXIS_VARNAMES)
+
+    ds["grid"] = xr.DataArray(
+        0,
+        attrs=sgrid.SGrid2DMetadata(
+            cf_role="grid_topology",
+            topology_dimension=2,
+            node_dimensions=("M", "N"),
+            node_coordinates=("lon", "lat"),
+            face_dimensions=(
+                sgrid.FaceNodePadding("X", "M", sgrid.Padding.LOW),
+                sgrid.FaceNodePadding("Y", "N", sgrid.Padding.LOW),
+            ),
+            vertical_dimensions=(sgrid.FaceNodePadding("Z", "LAYER", sgrid.Padding.HIGH),),
+        ).to_attrs(),
+    )
+
+    return ds
+
+
 # Known vertical dimension mappings by model
 _FESOM2_VERTICAL_DIMS = {"interface": "nz", "center": "nz1"}
 _ICON_VERTICAL_DIMS = {"interface": "depth_2", "center": "depth"}

tests/test_convert.py

@@ -159,6 +159,16 @@ def test_convert_copernicusmarine_no_currents(caplog):
     assert caplog.text == ""
 
 
+def test_convert_structured_delft3d():
+    ds = open_remote_dataset("Delft3D_data/Rotterdam_tiny")
+    coords = ds[["XZETA", "YZETA", "SIGMA_C"]]
+    ds_fset = convert.delft3d_to_sgrid(fields={"U": ds["VELU"], "V": ds["VELV"]}, coords=coords)
+    fieldset = FieldSet.from_sgrid_conventions(ds_fset)
+    assert "U" in fieldset.fields
+    assert "V" in fieldset.fields
+    assert "UV" in fieldset.fields
+
+
 @pytest.mark.parametrize("ds", _COPERNICUS_DATASETS)
 def test_convert_copernicusmarine_no_logs(ds, caplog):
     ds = ds.copy()