diff --git a/docs/user_guide/examples/tutorial_delft3d.ipynb b/docs/user_guide/examples/tutorial_delft3d.ipynb index 18f9c9d5b..a83329f92 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" @@ -40,13 +41,263 @@ "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." + "## 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, 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", + "\n", + "We start by opening the raw file." ] }, { - "cell_type": "markdown", + "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 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", + "node_y = raw_ds[\"Mesh_node_y\"].values\n", + "# 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", + "]" + ] + }, + { + "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 = [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", + " 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. First, we split each grid cell into triangles using a fan triangulation, keeping track of which original grid cell each triangle came from." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "tri_faces = []\n", + "tri_parent_face = []\n", + "\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)" + ] + }, + { + "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", + " 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", + ")" + ] + }, + { + "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": "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "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": "15", + "metadata": {}, + "source": [ + "We can then run a simple Parcels simulation with the regridded version of the model output." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "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": "17", + "metadata": {}, + "source": [ + "Finally, we plot the trajectories of the particles on top of the flow mesh." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "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", + "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": "19", + "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": "20", + "metadata": {}, + "source": [ + "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": "21", + "metadata": {}, "source": [ "```{note}\n", "While Delft3D provides velocities on a CGrid, we currently can't use the {py:obj}`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 {py:obj}`parcels.interpolators.XFreeslip` interpolator instead.\n", @@ -56,7 +307,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -76,7 +327,7 @@ }, { "cell_type": "markdown", - "id": "6", + "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." @@ -85,7 +336,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -96,7 +347,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", @@ -123,7 +374,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "25", "metadata": {}, "outputs": [], "source": [ @@ -168,7 +419,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Parcels:docs (3.14.6)", + "display_name": "Parcels:default (3.14.8)", "language": "python", "name": "python3" }, @@ -182,7 +433,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.6" + "version": "3.14.8" }, "mystnb": { "execution_in_temp": true 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)),