From 993cdda93a7634614baee530fc5cd7a0114181f5 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 17 Dec 2025 12:07:51 +0100 Subject: [PATCH 01/12] Creating a tutorial_nestedgrids Using idealised flow fields --- docs/examples/tutorial_nestedgrids.ipynb | 419 +++++++++++++++++++++++ docs/user_guide/index.md | 1 + pixi.toml | 1 + 3 files changed, 421 insertions(+) create mode 100644 docs/examples/tutorial_nestedgrids.ipynb diff --git a/docs/examples/tutorial_nestedgrids.ipynb b/docs/examples/tutorial_nestedgrids.ipynb new file mode 100644 index 000000000..3a5b05179 --- /dev/null +++ b/docs/examples/tutorial_nestedgrids.ipynb @@ -0,0 +1,419 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# 🖥️ Nested Grids\n", + "\n", + "In some applications, you may have access to fields on different grids that each cover only part of the region of interest. Then, you would like to combine them all together. You may also have a grid covering the entire region and another one only covering part of it, but with a higher resolution. The set of those grids form what we call nested grids.\n", + "\n", + "In Parcels v4, we can use the new `uxarray` integration to determine in which grid a particle is located. We will demonstrate how to set up a simulation with multiple nested grids, and how to handle particle transitions between these grids.\n", + "\n", + "\n", + "This tutorial shows how to use these Nested Fields with a very idealised example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import matplotlib.tri as mtri\n", + "import numpy as np\n", + "import uxarray as ux\n", + "import xarray as xr\n", + "from shapely.geometry import MultiPoint, Point, Polygon\n", + "from triangle import triangulate\n", + "\n", + "import parcels" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Setting up the individual Grids and Fields\n", + "\n", + "First define a zonal and meridional velocity field defined on a high resolution (dx = 100m) 2kmx2km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 * cos(lon / 200 * pi / 2) m/s." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "from parcels._datasets.structured.generated import simple_UV_dataset\n", + "\n", + "xdim = 10\n", + "ydim = 10\n", + "\n", + "ds_in = [None] * 2\n", + "grid_polygons = [None] * 2\n", + "\n", + "grid_polygons[0] = Polygon([(0, 0), (2000, 0), (2000, 2000), (0, 2000)])\n", + "ds_in[0] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"flat\").isel(time=0, depth=0)\n", + "\n", + "ds_in[0][\"lon\"][:] = np.linspace(\n", + " grid_polygons[0].bounds[0], grid_polygons[0].bounds[2], xdim\n", + ")\n", + "ds_in[0][\"lat\"][:] = np.linspace(\n", + " grid_polygons[0].bounds[1], grid_polygons[0].bounds[3], ydim\n", + ")\n", + "lon_g, lat_g = np.meshgrid(ds_in[0][\"lon\"], ds_in[0][\"lat\"])\n", + "ds_in[0][\"U\"][:] = 1.0\n", + "ds_in[0][\"V\"][:] = np.cos(lon_g / 200 * np.pi / 2)" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "Now define another velocity field on a low resolution (dx = 2km) 20kmx4km grid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "xdim = 40\n", + "ydim = 20\n", + "\n", + "grid_polygons[1] = Polygon(\n", + " [(-2000, -1000), (18000, -1000), (18000, 3000), (-2000, 3000)]\n", + ")\n", + "ds_in[1] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"flat\").isel(time=0, depth=0)\n", + "\n", + "ds_in[1][\"lon\"][:] = np.linspace(\n", + " grid_polygons[1].bounds[0], grid_polygons[1].bounds[2], xdim\n", + ")\n", + "ds_in[1][\"lat\"][:] = lat = np.linspace(\n", + " grid_polygons[1].bounds[1], grid_polygons[1].bounds[3], ydim\n", + ")\n", + "lon_g, lat_g = np.meshgrid(ds_in[1][\"lon\"], ds_in[1][\"lat\"])\n", + "ds_in[1][\"U\"][:] = 1.0\n", + "ds_in[1][\"V\"][:] = np.cos(lon_g / 200 * np.pi / 2)" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "Plot the two velocity fields on top of each other, indicating the grid boundaries in red." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "n_grids = len(ds_in)\n", + "fig, ax = plt.subplots()\n", + "\n", + "for i in range(n_grids):\n", + " ds = ds_in[i]\n", + " ax.quiver(ds.lon, ds.lat, ds.U, ds.V)\n", + " ax.plot(\n", + " np.append(grid_polygons[i].exterior.xy[0], grid_polygons[i].exterior.xy[0][0]),\n", + " np.append(grid_polygons[i].exterior.xy[1], grid_polygons[i].exterior.xy[1][0]),\n", + " \"-r\",\n", + " lw=2,\n", + " )\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Creating a Delaunay triangulation of the nested Grids\n", + "\n", + "Now comes the important part: we need to create a Delaunay triangulation of the nested Grids, so that we can efficiently determine in which Grid a particle is located at any given time. We use the `triangle` package to perform the triangulation, and `shapely` to handle the geometric operations. \n", + "\n", + "Note that we need to keep the edges of the polygons in the triangulation, so we need a [constrained (PSLG) Delaunay triangulation](https://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation).\n", + "\n", + "The result is a set of triangles covering the nested Grids, which we can use to determine in which Grid a particle is located at any given time. It is important that the list of polygons is ordered from smallest to largest Grid, so that triangles in overlapping areas are assigned to the correct Grid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from shapely.geometry import Point, Polygon\n", + "from triangle import triangulate\n", + "\n", + "\n", + "def constrained_triangulate_keep_edges(shapely_polys):\n", + " \"\"\"Constrained triangulation while keeping polygon edges.\n", + "\n", + " Accepts a list of shapely.geometry.Polygon objects.\n", + "\n", + " Returns:\n", + " pts (P,2) array of vertex coordinates,\n", + " tris (M,3) array of triangle vertex indices,\n", + " face_poly (M,) mapping each triangle to the polygon index or -1.\n", + " \"\"\"\n", + " # normalize input: convert shapely Polygons to Nx2 numpy arrays\n", + " np_polys = []\n", + " for p in shapely_polys:\n", + " coords = np.asarray(p.exterior.coords)\n", + " if coords.shape[0] > 1 and np.allclose(coords[0], coords[-1]):\n", + " coords = coords[:-1]\n", + " np_polys.append(coords.astype(float))\n", + "\n", + " # build vertices + segments for Triangle PSLG\n", + " verts = []\n", + " segments = []\n", + " offset = 0\n", + " for poly in np_polys:\n", + " Ni = len(poly)\n", + " verts.extend(poly.tolist())\n", + " segments.extend([[offset + j, offset + ((j + 1) % Ni)] for j in range(Ni)])\n", + " offset += Ni\n", + " verts = np.asarray(verts, dtype=float)\n", + " segments = np.asarray(segments, dtype=int)\n", + "\n", + " mode = \"p\" # \"p\" = PSLG (constrained triangulation)\n", + " B = triangulate({\"vertices\": verts, \"segments\": segments}, mode)\n", + "\n", + " pts = B[\"vertices\"]\n", + " tris = B[\"triangles\"].astype(int)\n", + "\n", + " # assign triangles to polygons using centroid test\n", + " centers = pts[tris].mean(axis=1)\n", + " face_poly = np.full(len(tris), -1, dtype=int)\n", + " for ti, c in enumerate(centers):\n", + " for ip in range(len(shapely_polys)):\n", + " if shapely_polys[ip].contains(Point(c)):\n", + " face_poly[ti] = ip\n", + " break\n", + "\n", + " return pts, tris, face_poly" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "We can then run the triangulation and plot the resulting triangles to verify that they correctly cover the nested Grids." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "points, face_tris, face_poly = constrained_triangulate_keep_edges(grid_polygons)\n", + "\n", + "fig, ax = plt.subplots()\n", + "for i in range(n_grids)[::-1]:\n", + " tris = face_tris[face_poly == i]\n", + " ax.triplot(points[:, 0], points[:, 1], tris, label=f\"Nest {i}\")\n", + "ax.scatter(points[:, 0], points[:, 1], s=10, c=\"k\")\n", + "ax.set_aspect(\"equal\")\n", + "ax.legend()\n", + "\n", + "# reverse legend labels\n", + "handles, labels = ax.get_legend_handles_labels()\n", + "ax.legend(handles[::-1], labels[::-1])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "Then, we convert the triangulation into a Parcels FieldSet using `Parcels.UxGrid()`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "# build an xarray dataset compatible with UGRID / uxarray\n", + "n_node = points.shape[0]\n", + "n_face = face_tris.shape[0]\n", + "n_max_face_nodes = face_tris.shape[1]\n", + "\n", + "ds_tri = xr.Dataset(\n", + " {\n", + " \"node_lon\": ((\"n_node\",), points[:, 0]),\n", + " \"node_lat\": ((\"n_node\",), points[:, 1]),\n", + " \"face_node_connectivity\": ((\"n_face\", \"n_max_face_nodes\"), face_tris),\n", + " \"face_polygon\": (\n", + " (\n", + " \"time\",\n", + " \"nz\",\n", + " \"n_face\",\n", + " ),\n", + " face_poly[np.newaxis, np.newaxis, :],\n", + " {\n", + " \"long_name\": \"Nested grid ID\",\n", + " \"location\": \"face\",\n", + " \"mesh\": \"delaunay\",\n", + " },\n", + " ),\n", + " },\n", + " coords={\n", + " \"time\": np.array([np.timedelta64(0, \"ns\")]),\n", + " \"nz\": np.array([0]),\n", + " \"n_node\": np.arange(n_node),\n", + " \"n_face\": np.arange(n_face),\n", + " },\n", + " attrs={\"Conventions\": \"UGRID-1.0\"},\n", + ")\n", + "\n", + "uxda = ux.UxDataArray(ds_tri[\"face_polygon\"], uxgrid=ux.Grid(ds_tri))\n", + "\n", + "NestID = parcels.Field(\n", + " \"NestID\",\n", + " uxda,\n", + " parcels.UxGrid(uxda.uxgrid, z=uxda[\"nz\"], mesh=\"flat\"),\n", + " interp_method=parcels.interpolators.UxPiecewiseConstantFace,\n", + ")\n", + "fieldset = parcels.FieldSet([NestID])" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "We can confirm that the FieldSet has been created correctly by running a Parcels simulation where particles sample the `NestID` field, which indicates in which nest each particle is located at any given time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "X, Y = np.meshgrid(\n", + " np.linspace(-1000, 5000, 10),\n", + " np.linspace(-500, 2500, 10),\n", + ")\n", + "\n", + "NestParticle = parcels.Particle.add_variable(parcels.Variable(\"nestID\", dtype=np.int32))\n", + "pset = parcels.ParticleSet(\n", + " fieldset, pclass=NestParticle, lon=X.flatten(), lat=Y.flatten()\n", + ")\n", + "\n", + "\n", + "def SampleNestID(particles, fieldset):\n", + " particles.nestID = fieldset.NestID[particles]\n", + "\n", + "\n", + "pset.execute(SampleNestID, runtime=np.timedelta64(1, \"s\"), dt=np.timedelta64(1, \"s\"))" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "Indeed, the visualisation below shows that particles correctly identify the nest they are in based on their location." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n", + "\n", + "triang = mtri.Triangulation(\n", + " uxda.uxgrid.node_lon.values,\n", + " uxda.uxgrid.node_lat.values,\n", + " triangles=uxda.uxgrid.face_node_connectivity.values,\n", + ")\n", + "\n", + "plot_args = {\n", + " \"cmap\": \"viridis\",\n", + " \"edgecolors\": \"k\",\n", + " \"linewidth\": 0.5,\n", + " \"vmin\": 0,\n", + " \"vmax\": 2.0,\n", + "}\n", + "ax.tripcolor(\n", + " triang, facecolors=np.squeeze(uxda[0, :].values), shading=\"flat\", **plot_args\n", + ")\n", + "ax.scatter(pset.lon, pset.lat, c=pset.nestID, **plot_args)\n", + "ax.set_aspect(\"equal\")\n", + "ax.set_title(\"Nesting visualisation (triangulation and interpolated particle values)\")\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "**Note for Joe:** Here it goes wrong. I think it's because of the UxGrid, wrapping around the antimeridian. See next cell." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "print(ds_tri.node_lon.values)\n", + "print(ux.Grid(ds_tri).node_lon.values)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "docs", + "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.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index a46289455..419594162 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -31,6 +31,7 @@ TODO: Add links to Reference API throughout examples/explanation_grids.md examples/tutorial_nemo_curvilinear.ipynb examples/tutorial_unitconverters.ipynb +examples/tutorial_nestedgrids.ipynb ``` diff --git a/pixi.toml b/pixi.toml index 7b6fafac1..3e28cc342 100644 --- a/pixi.toml +++ b/pixi.toml @@ -30,6 +30,7 @@ xgcm = ">=0.9.0" cf_xarray = ">=0.8.6" cftime = ">=1.6.3" pooch = ">=1.8.0" +py-triangle = ">=20250106,<20250107" [dependencies] parcels = { path = "." } From 8b98f1f63af363b6a121d11b9e3b5f5ec24eb2e9 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 17 Dec 2025 12:42:09 +0100 Subject: [PATCH 02/12] move tutorial to right location --- docs/{ => user_guide}/examples/tutorial_nestedgrids.ipynb | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/{ => user_guide}/examples/tutorial_nestedgrids.ipynb (100%) diff --git a/docs/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb similarity index 100% rename from docs/examples/tutorial_nestedgrids.ipynb rename to docs/user_guide/examples/tutorial_nestedgrids.ipynb From 37d6cf329030172ef0e26df451c567a4312f9410 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 17 Dec 2025 16:24:09 +0100 Subject: [PATCH 03/12] Using spherical mesh for nestedgrid tutorial --- .../examples/tutorial_nestedgrids.ipynb | 81 +++++++------------ 1 file changed, 29 insertions(+), 52 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index 3a5b05179..efecae888 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -30,7 +30,8 @@ "from shapely.geometry import MultiPoint, Point, Polygon\n", "from triangle import triangulate\n", "\n", - "import parcels" + "import parcels\n", + "from parcels._datasets.structured.generated import simple_UV_dataset" ] }, { @@ -40,7 +41,7 @@ "source": [ "## Setting up the individual Grids and Fields\n", "\n", - "First define a zonal and meridional velocity field defined on a high resolution (dx = 100m) 2kmx2km grid with a flat mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 * cos(lon / 200 * pi / 2) m/s." + "First define a zonal and meridional velocity field defined on a high resolution mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 * cos(lon * pi / 2) m/s." ] }, { @@ -50,26 +51,26 @@ "metadata": {}, "outputs": [], "source": [ - "from parcels._datasets.structured.generated import simple_UV_dataset\n", - "\n", "xdim = 10\n", "ydim = 10\n", "\n", "ds_in = [None] * 2\n", "grid_polygons = [None] * 2\n", "\n", - "grid_polygons[0] = Polygon([(0, 0), (2000, 0), (2000, 2000), (0, 2000)])\n", - "ds_in[0] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"flat\").isel(time=0, depth=0)\n", + "grid_polygons[0] = np.array([(0, 0), (10, 0), (10, 10), (0, 10)])\n", + "ds_in[0] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"spherical\").isel(\n", + " time=0, depth=0\n", + ")\n", "\n", "ds_in[0][\"lon\"][:] = np.linspace(\n", - " grid_polygons[0].bounds[0], grid_polygons[0].bounds[2], xdim\n", + " grid_polygons[0][:, 0].min(), grid_polygons[0][:, 0].max(), xdim\n", ")\n", "ds_in[0][\"lat\"][:] = np.linspace(\n", - " grid_polygons[0].bounds[1], grid_polygons[0].bounds[3], ydim\n", + " grid_polygons[0][:, 1].min(), grid_polygons[0][:, 1].max(), ydim\n", ")\n", "lon_g, lat_g = np.meshgrid(ds_in[0][\"lon\"], ds_in[0][\"lat\"])\n", "ds_in[0][\"U\"][:] = 1.0\n", - "ds_in[0][\"V\"][:] = np.cos(lon_g / 200 * np.pi / 2)" + "ds_in[0][\"V\"][:] = np.cos(lon_g * np.pi / 2)" ] }, { @@ -77,7 +78,7 @@ "id": "4", "metadata": {}, "source": [ - "Now define another velocity field on a low resolution (dx = 2km) 20kmx4km grid." + "Now define another velocity field on a lower resolution grid." ] }, { @@ -90,20 +91,20 @@ "xdim = 40\n", "ydim = 20\n", "\n", - "grid_polygons[1] = Polygon(\n", - " [(-2000, -1000), (18000, -1000), (18000, 3000), (-2000, 3000)]\n", + "grid_polygons[1] = np.array([(-10, -20), (60, -20), (60, 40), (-10, 40)])\n", + "ds_in[1] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"spherical\").isel(\n", + " time=0, depth=0\n", ")\n", - "ds_in[1] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"flat\").isel(time=0, depth=0)\n", "\n", "ds_in[1][\"lon\"][:] = np.linspace(\n", - " grid_polygons[1].bounds[0], grid_polygons[1].bounds[2], xdim\n", + " grid_polygons[1][:, 0].min(), grid_polygons[1][:, 0].max(), xdim\n", ")\n", - "ds_in[1][\"lat\"][:] = lat = np.linspace(\n", - " grid_polygons[1].bounds[1], grid_polygons[1].bounds[3], ydim\n", + "ds_in[1][\"lat\"][:] = np.linspace(\n", + " grid_polygons[1][:, 1].min(), grid_polygons[1][:, 1].max(), ydim\n", ")\n", "lon_g, lat_g = np.meshgrid(ds_in[1][\"lon\"], ds_in[1][\"lat\"])\n", "ds_in[1][\"U\"][:] = 1.0\n", - "ds_in[1][\"V\"][:] = np.cos(lon_g / 200 * np.pi / 2)" + "ds_in[1][\"V\"][:] = np.cos(lon_g * 5 * np.pi / 2)" ] }, { @@ -127,9 +128,10 @@ "for i in range(n_grids):\n", " ds = ds_in[i]\n", " ax.quiver(ds.lon, ds.lat, ds.U, ds.V)\n", + " poly = grid_polygons[i]\n", " ax.plot(\n", - " np.append(grid_polygons[i].exterior.xy[0], grid_polygons[i].exterior.xy[0][0]),\n", - " np.append(grid_polygons[i].exterior.xy[1], grid_polygons[i].exterior.xy[1][0]),\n", + " np.append(poly[:, 0], poly[0, 0]),\n", + " np.append(poly[:, 1], poly[1, 1]),\n", " \"-r\",\n", " lw=2,\n", " )\n", @@ -164,29 +166,22 @@ "from triangle import triangulate\n", "\n", "\n", - "def constrained_triangulate_keep_edges(shapely_polys):\n", + "def constrained_triangulate_keep_edges(polygons):\n", " \"\"\"Constrained triangulation while keeping polygon edges.\n", "\n", - " Accepts a list of shapely.geometry.Polygon objects.\n", + " Args:\n", + " polygons: list of (Ni,2) numpy arrays (polygon boundary points in order)\n", "\n", " Returns:\n", " pts (P,2) array of vertex coordinates,\n", " tris (M,3) array of triangle vertex indices,\n", " face_poly (M,) mapping each triangle to the polygon index or -1.\n", " \"\"\"\n", - " # normalize input: convert shapely Polygons to Nx2 numpy arrays\n", - " np_polys = []\n", - " for p in shapely_polys:\n", - " coords = np.asarray(p.exterior.coords)\n", - " if coords.shape[0] > 1 and np.allclose(coords[0], coords[-1]):\n", - " coords = coords[:-1]\n", - " np_polys.append(coords.astype(float))\n", - "\n", " # build vertices + segments for Triangle PSLG\n", " verts = []\n", " segments = []\n", " offset = 0\n", - " for poly in np_polys:\n", + " for poly in polygons:\n", " Ni = len(poly)\n", " verts.extend(poly.tolist())\n", " segments.extend([[offset + j, offset + ((j + 1) % Ni)] for j in range(Ni)])\n", @@ -201,6 +196,7 @@ " tris = B[\"triangles\"].astype(int)\n", "\n", " # assign triangles to polygons using centroid test\n", + " shapely_polys = [Polygon(p) for p in polygons]\n", " centers = pts[tris].mean(axis=1)\n", " face_poly = np.full(len(tris), -1, dtype=int)\n", " for ti, c in enumerate(centers):\n", @@ -296,7 +292,7 @@ "NestID = parcels.Field(\n", " \"NestID\",\n", " uxda,\n", - " parcels.UxGrid(uxda.uxgrid, z=uxda[\"nz\"], mesh=\"flat\"),\n", + " parcels.UxGrid(uxda.uxgrid, z=uxda[\"nz\"], mesh=\"spherical\"),\n", " interp_method=parcels.interpolators.UxPiecewiseConstantFace,\n", ")\n", "fieldset = parcels.FieldSet([NestID])" @@ -318,8 +314,8 @@ "outputs": [], "source": [ "X, Y = np.meshgrid(\n", - " np.linspace(-1000, 5000, 10),\n", - " np.linspace(-500, 2500, 10),\n", + " np.linspace(1, 9, 5),\n", + " np.linspace(1, 9, 5),\n", ")\n", "\n", "NestParticle = parcels.Particle.add_variable(parcels.Variable(\"nestID\", dtype=np.int32))\n", @@ -374,25 +370,6 @@ "plt.tight_layout()\n", "plt.show()" ] - }, - { - "cell_type": "markdown", - "id": "18", - "metadata": {}, - "source": [ - "**Note for Joe:** Here it goes wrong. I think it's because of the UxGrid, wrapping around the antimeridian. See next cell." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "19", - "metadata": {}, - "outputs": [], - "source": [ - "print(ds_tri.node_lon.values)\n", - "print(ux.Grid(ds_tri).node_lon.values)" - ] } ], "metadata": { From 41251561b3bbf173e5c1defaa954e3557b0f17cc Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 17 Dec 2025 17:09:56 +0100 Subject: [PATCH 04/12] Adding advection code to nestedgrids tutorial --- .../examples/tutorial_nestedgrids.ipynb | 159 +++++++++++++++++- 1 file changed, 153 insertions(+), 6 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index efecae888..d3c892a65 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -41,7 +41,7 @@ "source": [ "## Setting up the individual Grids and Fields\n", "\n", - "First define a zonal and meridional velocity field defined on a high resolution mesh. The zonal velocity is uniform and 1 m/s, and the meridional velocity is equal to 0.5 * cos(lon * pi / 2) m/s." + "First define a zonal and meridional velocity field defined on a high resolution mesh. Both the zonal and meridional velocity are uniform, with velocities of 1 m/s and 0 m/s respectively." ] }, { @@ -70,7 +70,7 @@ ")\n", "lon_g, lat_g = np.meshgrid(ds_in[0][\"lon\"], ds_in[0][\"lat\"])\n", "ds_in[0][\"U\"][:] = 1.0\n", - "ds_in[0][\"V\"][:] = np.cos(lon_g * np.pi / 2)" + "ds_in[0][\"V\"][:] = 0.0" ] }, { @@ -78,7 +78,7 @@ "id": "4", "metadata": {}, "source": [ - "Now define another velocity field on a lower resolution grid." + "Now define another velocity field on a coarser resolution grid. The zonal velocity is the same as for the finer resolution grid, but the meridional velocity is now a cosine as a function of longitude." ] }, { @@ -104,7 +104,7 @@ ")\n", "lon_g, lat_g = np.meshgrid(ds_in[1][\"lon\"], ds_in[1][\"lat\"])\n", "ds_in[1][\"U\"][:] = 1.0\n", - "ds_in[1][\"V\"][:] = np.cos(lon_g * 5 * np.pi / 2)" + "ds_in[1][\"V\"][:] = np.cos(lon_g / 5 * np.pi / 2)" ] }, { @@ -292,7 +292,8 @@ "NestID = parcels.Field(\n", " \"NestID\",\n", " uxda,\n", - " parcels.UxGrid(uxda.uxgrid, z=uxda[\"nz\"], mesh=\"spherical\"),\n", + " # TODO note that here we need to use mesh=\"flat\" otherwise the hashing doesn't work. See https://github.com/Parcels-code/Parcels/pull/2439#discussion_r2627664010\n", + " parcels.UxGrid(uxda.uxgrid, z=uxda[\"nz\"], mesh=\"flat\"),\n", " interp_method=parcels.interpolators.UxPiecewiseConstantFace,\n", ")\n", "fieldset = parcels.FieldSet([NestID])" @@ -347,6 +348,7 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n", + "cmap = plt.get_cmap(\"viridis\", n_grids)\n", "\n", "triang = mtri.Triangulation(\n", " uxda.uxgrid.node_lon.values,\n", @@ -355,7 +357,7 @@ ")\n", "\n", "plot_args = {\n", - " \"cmap\": \"viridis\",\n", + " \"cmap\": cmap,\n", " \"edgecolors\": \"k\",\n", " \"linewidth\": 0.5,\n", " \"vmin\": 0,\n", @@ -370,6 +372,151 @@ "plt.tight_layout()\n", "plt.show()" ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "## Advecting particles with nest transitions\n", + "\n", + "We can now set up a particle advection simulation using the nested grids. We first combine all the Fields into a single FieldSet. \n", + "\n", + "We rename the individual Fields by appending the nest index to their names, so that we can easily identify which Field belongs to which nest. We also add the `NestID` Field to the FieldSet (note that Parcels v4 supports combining structured and unstructured Fields into one FieldSet, which is very convenient for this usecase)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "fields = [NestID]\n", + "for i, ds in enumerate(ds_in):\n", + " # TODO : use FieldSet.from_sgrid_convetion here once #2437 is merged\n", + " grid = parcels.XGrid.from_dataset(ds, mesh=\"flat\")\n", + " U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", + " V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", + " UV = parcels.VectorField(\"UV\", U, V)\n", + "\n", + " fset = parcels.FieldSet([U, V, UV])\n", + "\n", + " # Make sure to attach the correct unit converter (since field-names have been changed)\n", + " fset.U.units = parcels.GeographicPolar()\n", + " fset.V.units = parcels.Geographic()\n", + "\n", + " #\n", + " for fld in fset.fields.values():\n", + " fld.name = f\"{fld.name}{i}\"\n", + " fields.append(fld)\n", + "fieldset = parcels.FieldSet(fields)" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "We then define a custom Advection kernel that advects particles using the appropriate velocity Field based on the `NestID` at the particle's location. Note that for simplicity, we use Eulerian advection here, but in a real-world application you would typically use a higher-order scheme." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "def AdvectEE_Nests(particles, fieldset):\n", + " particles.nestID = fieldset.NestID[particles]\n", + "\n", + " # TODO because of KernelParticle bug (GH #2143), we need to copy lon/lat/time to local variables\n", + " time = particles.time\n", + " z = particles.z\n", + " lat = particles.lat\n", + " lon = particles.lon\n", + " u = np.zeros_like(particles.lon)\n", + " v = np.zeros_like(particles.lat)\n", + "\n", + " unique_ids = np.unique(particles.nestID)\n", + " for nid in unique_ids:\n", + " mask = particles.nestID == nid\n", + " UVField = getattr(fieldset, f\"UV{nid}\")\n", + " (u[mask], v[mask]) = UVField[time[mask], z[mask], lat[mask], lon[mask]]\n", + "\n", + " particles.dlon += u * particles.dt\n", + " particles.dlat += v * particles.dt\n", + "\n", + " # TODO particle states have to be updated manually because UVField is not called with `particles` argument (becaise of GH #2143)\n", + " particles.state = np.where(\n", + " np.isnan(u) | np.isnan(v),\n", + " parcels.StatusCode.ErrorInterpolation,\n", + " particles.state,\n", + " )\n", + "\n", + "\n", + "def DeleteErrorParticles(particles, fieldset):\n", + " any_error = particles.state >= 50 # This captures all Errors\n", + " particles[any_error].state = parcels.StatusCode.Delete\n", + "\n", + "\n", + "lat = np.linspace(-15, 25, 10)\n", + "lon = np.full(len(lat), -5)\n", + "\n", + "pset = parcels.ParticleSet(fieldset, pclass=NestParticle, lon=lon, lat=lat)\n", + "ofile = parcels.ParticleFile(\"nest_particles.zarr\", outputdt=np.timedelta64(1, \"D\"))\n", + "pset.execute(\n", + " [AdvectEE_Nests, DeleteErrorParticles],\n", + " runtime=np.timedelta64(40, \"D\"),\n", + " dt=np.timedelta64(1, \"D\"),\n", + " output_file=ofile,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n", + "\n", + "ds_out = xr.open_zarr(\"nest_particles.zarr\")\n", + "\n", + "plt.plot(ds_out.lon.T, ds_out.lat.T, \"k\", linewidth=0.5)\n", + "sc = ax.scatter(ds_out.lon, ds_out.lat, c=ds_out.nestID, s=4, cmap=cmap, vmin=0, vmax=2)\n", + "xl, yl = ax.get_xlim(), ax.get_ylim()\n", + "\n", + "for i in range(n_grids):\n", + " poly = grid_polygons[i]\n", + " ax.plot(\n", + " np.append(poly[:, 0], poly[0, 0]),\n", + " np.append(poly[:, 1], poly[1, 1]),\n", + " \"-r\",\n", + " lw=2,\n", + " )\n", + "ax.set_xlim(xl)\n", + "ax.set_ylim(yl)\n", + "ax.set_aspect(\"equal\")\n", + "\n", + "cbar = plt.colorbar(sc, ticks=[0, 1, 2], ax=ax)\n", + "cbar.set_label(\"Nest ID\")\n", + "ax.set_title(\"Particle advection through nests\")\n", + "plt.tight_layout\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "Indeed, we can see that the particles follow the cosine oscillation pattern in the coarser grid, and move purely zonally in the finer grid.\n", + "\n", + "Note that in the plot above the particles at higher latitudes move farther in longitude-space (at a constant zonal velocity) because of the curvature of the earth, so that is expected." + ] } ], "metadata": { From 7169a3a6503a44b30e335a9e558545138e00ed85 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 17 Dec 2025 17:49:26 +0100 Subject: [PATCH 05/12] Adding a third grid to nests --- .../examples/tutorial_nestedgrids.ipynb | 176 +++++++++++------- 1 file changed, 104 insertions(+), 72 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index d3c892a65..b4c49d79f 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -27,6 +27,7 @@ "import numpy as np\n", "import uxarray as ux\n", "import xarray as xr\n", + "from matplotlib.colors import ListedColormap\n", "from shapely.geometry import MultiPoint, Point, Polygon\n", "from triangle import triangulate\n", "\n", @@ -39,9 +40,7 @@ "id": "2", "metadata": {}, "source": [ - "## Setting up the individual Grids and Fields\n", - "\n", - "First define a zonal and meridional velocity field defined on a high resolution mesh. Both the zonal and meridional velocity are uniform, with velocities of 1 m/s and 0 m/s respectively." + "## Setting up the individual Grids and Fields" ] }, { @@ -51,26 +50,19 @@ "metadata": {}, "outputs": [], "source": [ - "xdim = 10\n", - "ydim = 10\n", + "ds_in = []\n", + "grid_polygons = []\n", "\n", - "ds_in = [None] * 2\n", - "grid_polygons = [None] * 2\n", "\n", - "grid_polygons[0] = np.array([(0, 0), (10, 0), (10, 10), (0, 10)])\n", - "ds_in[0] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"spherical\").isel(\n", - " time=0, depth=0\n", - ")\n", + "def setup_nestedFields(polygon):\n", + " xdim, ydim = 20, 20\n", "\n", - "ds_in[0][\"lon\"][:] = np.linspace(\n", - " grid_polygons[0][:, 0].min(), grid_polygons[0][:, 0].max(), xdim\n", - ")\n", - "ds_in[0][\"lat\"][:] = np.linspace(\n", - " grid_polygons[0][:, 1].min(), grid_polygons[0][:, 1].max(), ydim\n", - ")\n", - "lon_g, lat_g = np.meshgrid(ds_in[0][\"lon\"], ds_in[0][\"lat\"])\n", - "ds_in[0][\"U\"][:] = 1.0\n", - "ds_in[0][\"V\"][:] = 0.0" + " ds = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"spherical\").isel(\n", + " time=0, depth=0\n", + " )\n", + " ds[\"lon\"][:] = np.linspace(polygon[:, 0].min(), polygon[:, 0].max(), xdim)\n", + " ds[\"lat\"][:] = np.linspace(polygon[:, 1].min(), polygon[:, 1].max(), ydim)\n", + " return ds" ] }, { @@ -78,7 +70,7 @@ "id": "4", "metadata": {}, "source": [ - "Now define another velocity field on a coarser resolution grid. The zonal velocity is the same as for the finer resolution grid, but the meridional velocity is now a cosine as a function of longitude." + "First define a zonal and meridional velocity field defined on a high resolution mesh. Both the zonal and meridional velocity are uniform, with velocities of 1 m/s and 0 m/s respectively." ] }, { @@ -88,23 +80,12 @@ "metadata": {}, "outputs": [], "source": [ - "xdim = 40\n", - "ydim = 20\n", + "polygon = np.array([(10, 15), (25, 10), (25, 25), (10, 35)])\n", + "grid_polygons.append(polygon)\n", "\n", - "grid_polygons[1] = np.array([(-10, -20), (60, -20), (60, 40), (-10, 40)])\n", - "ds_in[1] = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"spherical\").isel(\n", - " time=0, depth=0\n", - ")\n", - "\n", - "ds_in[1][\"lon\"][:] = np.linspace(\n", - " grid_polygons[1][:, 0].min(), grid_polygons[1][:, 0].max(), xdim\n", - ")\n", - "ds_in[1][\"lat\"][:] = np.linspace(\n", - " grid_polygons[1][:, 1].min(), grid_polygons[1][:, 1].max(), ydim\n", - ")\n", - "lon_g, lat_g = np.meshgrid(ds_in[1][\"lon\"], ds_in[1][\"lat\"])\n", - "ds_in[1][\"U\"][:] = 1.0\n", - "ds_in[1][\"V\"][:] = np.cos(lon_g / 5 * np.pi / 2)" + "ds_in.append(setup_nestedFields(polygon))\n", + "ds_in[-1][\"U\"][:] = 1.0\n", + "ds_in[-1][\"V\"][:] = 0.0" ] }, { @@ -112,7 +93,7 @@ "id": "6", "metadata": {}, "source": [ - "Plot the two velocity fields on top of each other, indicating the grid boundaries in red." + "Then define another set of zonal and meridional velocities on a slightly coarser grid. In this case, the meridional velocity is also 1 m/s." ] }, { @@ -121,19 +102,68 @@ "id": "7", "metadata": {}, "outputs": [], + "source": [ + "polygon = np.array([(0, -5), (35, 0), (35, 25), (0, 20)])\n", + "grid_polygons.append(polygon)\n", + "\n", + "ds_in.append(setup_nestedFields(polygon))\n", + "ds_in[-1][\"U\"][:] = 1.0\n", + "ds_in[-1][\"V\"][:] = 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "Now define another velocity field on a coarser resolution grid. The zonal velocity is the same as for the finer resolution grid, but the meridional velocity is now a cosine as a function of longitude." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "polygon = np.array([(-10, -20), (60, -20), (60, 40), (-10, 40)])\n", + "grid_polygons.append(polygon)\n", + "\n", + "ds_in.append(setup_nestedFields(polygon))\n", + "lon_g, lat_g = np.meshgrid(ds_in[-1][\"lon\"], ds_in[-1][\"lat\"])\n", + "ds_in[-1][\"U\"][:] = 1.0\n", + "ds_in[-1][\"V\"][:] = np.cos(lon_g / 5 * np.pi / 2)" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "Plot the velocity fields on top of each other, indicating the grid boundaries in red." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], "source": [ "n_grids = len(ds_in)\n", + "cmap = ListedColormap([f\"C{i}\" for i in range(n_grids)])\n", + "\n", "fig, ax = plt.subplots()\n", "\n", "for i in range(n_grids):\n", " ds = ds_in[i]\n", - " ax.quiver(ds.lon, ds.lat, ds.U, ds.V)\n", + " ax.quiver(ds.lon, ds.lat, ds.U, ds.V, color=f\"C{i}\")\n", " poly = grid_polygons[i]\n", " ax.plot(\n", " np.append(poly[:, 0], poly[0, 0]),\n", - " np.append(poly[:, 1], poly[1, 1]),\n", + " np.append(poly[:, 1], poly[0, 1]),\n", " \"-r\",\n", - " lw=2,\n", + " lw=1,\n", " )\n", "\n", "plt.tight_layout()\n", @@ -142,7 +172,15 @@ }, { "cell_type": "markdown", - "id": "8", + "id": "12", + "metadata": {}, + "source": [ + "Note in the plot above that the dataset domains are rectangular, but the polygons that will later define the Nested grid boundaries aren't necessarily. So we can even use this method to subset parts of a Grid." + ] + }, + { + "cell_type": "markdown", + "id": "13", "metadata": {}, "source": [ "## Creating a Delaunay triangulation of the nested Grids\n", @@ -157,7 +195,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -210,7 +248,7 @@ }, { "cell_type": "markdown", - "id": "10", + "id": "15", "metadata": {}, "source": [ "We can then run the triangulation and plot the resulting triangles to verify that they correctly cover the nested Grids." @@ -219,7 +257,7 @@ { "cell_type": "code", "execution_count": null, - "id": "11", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -228,7 +266,7 @@ "fig, ax = plt.subplots()\n", "for i in range(n_grids)[::-1]:\n", " tris = face_tris[face_poly == i]\n", - " ax.triplot(points[:, 0], points[:, 1], tris, label=f\"Nest {i}\")\n", + " ax.triplot(points[:, 0], points[:, 1], tris, label=f\"Nest {i}\", color=f\"C{i}\")\n", "ax.scatter(points[:, 0], points[:, 1], s=10, c=\"k\")\n", "ax.set_aspect(\"equal\")\n", "ax.legend()\n", @@ -241,7 +279,7 @@ }, { "cell_type": "markdown", - "id": "12", + "id": "17", "metadata": {}, "source": [ "Then, we convert the triangulation into a Parcels FieldSet using `Parcels.UxGrid()`." @@ -250,7 +288,7 @@ { "cell_type": "code", "execution_count": null, - "id": "13", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -301,7 +339,7 @@ }, { "cell_type": "markdown", - "id": "14", + "id": "19", "metadata": {}, "source": [ "We can confirm that the FieldSet has been created correctly by running a Parcels simulation where particles sample the `NestID` field, which indicates in which nest each particle is located at any given time." @@ -310,13 +348,13 @@ { "cell_type": "code", "execution_count": null, - "id": "15", + "id": "20", "metadata": {}, "outputs": [], "source": [ "X, Y = np.meshgrid(\n", - " np.linspace(1, 9, 5),\n", - " np.linspace(1, 9, 5),\n", + " np.linspace(-8, 58, 25),\n", + " np.linspace(-18, 38, 20),\n", ")\n", "\n", "NestParticle = parcels.Particle.add_variable(parcels.Variable(\"nestID\", dtype=np.int32))\n", @@ -334,7 +372,7 @@ }, { "cell_type": "markdown", - "id": "16", + "id": "21", "metadata": {}, "source": [ "Indeed, the visualisation below shows that particles correctly identify the nest they are in based on their location." @@ -343,12 +381,11 @@ { "cell_type": "code", "execution_count": null, - "id": "17", + "id": "22", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n", - "cmap = plt.get_cmap(\"viridis\", n_grids)\n", "\n", "triang = mtri.Triangulation(\n", " uxda.uxgrid.node_lon.values,\n", @@ -375,7 +412,7 @@ }, { "cell_type": "markdown", - "id": "18", + "id": "23", "metadata": {}, "source": [ "## Advecting particles with nest transitions\n", @@ -388,25 +425,20 @@ { "cell_type": "code", "execution_count": null, - "id": "19", + "id": "24", "metadata": {}, "outputs": [], "source": [ "fields = [NestID]\n", "for i, ds in enumerate(ds_in):\n", " # TODO : use FieldSet.from_sgrid_convetion here once #2437 is merged\n", - " grid = parcels.XGrid.from_dataset(ds, mesh=\"flat\")\n", + " grid = parcels.XGrid.from_dataset(ds, mesh=\"spherical\")\n", " U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", " V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", " UV = parcels.VectorField(\"UV\", U, V)\n", "\n", " fset = parcels.FieldSet([U, V, UV])\n", "\n", - " # Make sure to attach the correct unit converter (since field-names have been changed)\n", - " fset.U.units = parcels.GeographicPolar()\n", - " fset.V.units = parcels.Geographic()\n", - "\n", - " #\n", " for fld in fset.fields.values():\n", " fld.name = f\"{fld.name}{i}\"\n", " fields.append(fld)\n", @@ -415,7 +447,7 @@ }, { "cell_type": "markdown", - "id": "20", + "id": "25", "metadata": {}, "source": [ "We then define a custom Advection kernel that advects particles using the appropriate velocity Field based on the `NestID` at the particle's location. Note that for simplicity, we use Eulerian advection here, but in a real-world application you would typically use a higher-order scheme." @@ -424,7 +456,7 @@ { "cell_type": "code", "execution_count": null, - "id": "21", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -461,14 +493,14 @@ " particles[any_error].state = parcels.StatusCode.Delete\n", "\n", "\n", - "lat = np.linspace(-15, 25, 10)\n", + "lat = np.linspace(-17, 35, 10)\n", "lon = np.full(len(lat), -5)\n", "\n", "pset = parcels.ParticleSet(fieldset, pclass=NestParticle, lon=lon, lat=lat)\n", "ofile = parcels.ParticleFile(\"nest_particles.zarr\", outputdt=np.timedelta64(1, \"D\"))\n", "pset.execute(\n", " [AdvectEE_Nests, DeleteErrorParticles],\n", - " runtime=np.timedelta64(40, \"D\"),\n", + " runtime=np.timedelta64(60, \"D\"),\n", " dt=np.timedelta64(1, \"D\"),\n", " output_file=ofile,\n", ")" @@ -477,7 +509,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22", + "id": "27", "metadata": {}, "outputs": [], "source": [ @@ -489,13 +521,13 @@ "sc = ax.scatter(ds_out.lon, ds_out.lat, c=ds_out.nestID, s=4, cmap=cmap, vmin=0, vmax=2)\n", "xl, yl = ax.get_xlim(), ax.get_ylim()\n", "\n", - "for i in range(n_grids):\n", + "for i in range(n_grids - 1):\n", " poly = grid_polygons[i]\n", " ax.plot(\n", " np.append(poly[:, 0], poly[0, 0]),\n", - " np.append(poly[:, 1], poly[1, 1]),\n", + " np.append(poly[:, 1], poly[0, 1]),\n", " \"-r\",\n", - " lw=2,\n", + " lw=1,\n", " )\n", "ax.set_xlim(xl)\n", "ax.set_ylim(yl)\n", @@ -510,7 +542,7 @@ }, { "cell_type": "markdown", - "id": "23", + "id": "28", "metadata": {}, "source": [ "Indeed, we can see that the particles follow the cosine oscillation pattern in the coarser grid, and move purely zonally in the finer grid.\n", From 400996489ee1d3156c83c5eadf69cdc14785197a Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 17 Dec 2025 17:52:46 +0100 Subject: [PATCH 06/12] Minor clean-ups --- .../examples/tutorial_nestedgrids.ipynb | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index b4c49d79f..dd6d1e751 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -40,7 +40,9 @@ "id": "2", "metadata": {}, "source": [ - "## Setting up the individual Grids and Fields" + "## Setting up the individual Grids and Fields\n", + "\n", + "Before we start, define a helper function to quickly set up a nested dataset" ] }, { @@ -54,12 +56,9 @@ "grid_polygons = []\n", "\n", "\n", - "def setup_nestedFields(polygon):\n", + "def setup_nested_ds(polygon):\n", " xdim, ydim = 20, 20\n", - "\n", - " ds = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"spherical\").isel(\n", - " time=0, depth=0\n", - " )\n", + " ds = simple_UV_dataset(dims=(1, 1, ydim, xdim), mesh=\"spherical\")\n", " ds[\"lon\"][:] = np.linspace(polygon[:, 0].min(), polygon[:, 0].max(), xdim)\n", " ds[\"lat\"][:] = np.linspace(polygon[:, 1].min(), polygon[:, 1].max(), ydim)\n", " return ds" @@ -83,7 +82,7 @@ "polygon = np.array([(10, 15), (25, 10), (25, 25), (10, 35)])\n", "grid_polygons.append(polygon)\n", "\n", - "ds_in.append(setup_nestedFields(polygon))\n", + "ds_in.append(setup_nested_ds(polygon))\n", "ds_in[-1][\"U\"][:] = 1.0\n", "ds_in[-1][\"V\"][:] = 0.0" ] @@ -106,7 +105,7 @@ "polygon = np.array([(0, -5), (35, 0), (35, 25), (0, 20)])\n", "grid_polygons.append(polygon)\n", "\n", - "ds_in.append(setup_nestedFields(polygon))\n", + "ds_in.append(setup_nested_ds(polygon))\n", "ds_in[-1][\"U\"][:] = 1.0\n", "ds_in[-1][\"V\"][:] = 1.0" ] @@ -129,7 +128,7 @@ "polygon = np.array([(-10, -20), (60, -20), (60, 40), (-10, 40)])\n", "grid_polygons.append(polygon)\n", "\n", - "ds_in.append(setup_nestedFields(polygon))\n", + "ds_in.append(setup_nested_ds(polygon))\n", "lon_g, lat_g = np.meshgrid(ds_in[-1][\"lon\"], ds_in[-1][\"lat\"])\n", "ds_in[-1][\"U\"][:] = 1.0\n", "ds_in[-1][\"V\"][:] = np.cos(lon_g / 5 * np.pi / 2)" @@ -156,7 +155,7 @@ "fig, ax = plt.subplots()\n", "\n", "for i in range(n_grids):\n", - " ds = ds_in[i]\n", + " ds = ds_in[i].isel(time=0, depth=0)\n", " ax.quiver(ds.lon, ds.lat, ds.U, ds.V, color=f\"C{i}\")\n", " poly = grid_polygons[i]\n", " ax.plot(\n", From ea52655d632aba692a03fb3e50c91789a6e4fa46 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Wed, 17 Dec 2025 17:57:19 +0100 Subject: [PATCH 07/12] Suppressing pset.execute() output --- docs/user_guide/examples/tutorial_nestedgrids.ipynb | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index dd6d1e751..b2fac3d93 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -502,6 +502,7 @@ " runtime=np.timedelta64(60, \"D\"),\n", " dt=np.timedelta64(1, \"D\"),\n", " output_file=ofile,\n", + " verbose_progress=False,\n", ")" ] }, From e2ddeb173274a5e17c9d3be2bd0a516bcb1c6084 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 18 Dec 2025 07:52:24 +0100 Subject: [PATCH 08/12] Cleaning up text and code on nestedgrid tutorial --- .../examples/tutorial_nestedgrids.ipynb | 125 +++++++++--------- 1 file changed, 62 insertions(+), 63 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index b2fac3d93..40eb8315a 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -42,7 +42,7 @@ "source": [ "## Setting up the individual Grids and Fields\n", "\n", - "Before we start, define a helper function to quickly set up a nested dataset" + "First we define a helper function to quickly set up a nested dataset" ] }, { @@ -69,7 +69,7 @@ "id": "4", "metadata": {}, "source": [ - "First define a zonal and meridional velocity field defined on a high resolution mesh. Both the zonal and meridional velocity are uniform, with velocities of 1 m/s and 0 m/s respectively." + "Now we create a zonal and meridional velocity field defined on a small grid. Both the zonal (1m/s) and meridional (0 m/s) velocity are uniform." ] }, { @@ -92,7 +92,7 @@ "id": "6", "metadata": {}, "source": [ - "Then define another set of zonal and meridional velocities on a slightly coarser grid. In this case, the meridional velocity is also 1 m/s." + "Then, we define another set of zonal and meridional velocities on a slightly larger grid. In this case, the meridional velocity is also 1 m/s." ] }, { @@ -115,7 +115,7 @@ "id": "8", "metadata": {}, "source": [ - "Now define another velocity field on a coarser resolution grid. The zonal velocity is the same as for the finer resolution grid, but the meridional velocity is now a cosine as a function of longitude." + "Finally, we define another set of velocities on an even larger grid. The zonal velocity is the same as for the smaller grids, but the meridional velocity is now a cosine as a function of longitude." ] }, { @@ -139,7 +139,7 @@ "id": "10", "metadata": {}, "source": [ - "Plot the velocity fields on top of each other, indicating the grid boundaries in red." + "We plot the velocity fields on top of each other, indicating the grid boundaries in red." ] }, { @@ -174,7 +174,7 @@ "id": "12", "metadata": {}, "source": [ - "Note in the plot above that the dataset domains are rectangular, but the polygons that will later define the Nested grid boundaries aren't necessarily. So we can even use this method to subset parts of a Grid." + "Note, as seen in the plot above, that the dataset domains in this case are rectangular, but the polygons that will later define the nested Grid boundaries don't have to be. So we can even use this method to subset parts of a Grid." ] }, { @@ -186,7 +186,7 @@ "\n", "Now comes the important part: we need to create a Delaunay triangulation of the nested Grids, so that we can efficiently determine in which Grid a particle is located at any given time. We use the `triangle` package to perform the triangulation, and `shapely` to handle the geometric operations. \n", "\n", - "Note that we need to keep the edges of the polygons in the triangulation, so we need a [constrained (PSLG) Delaunay triangulation](https://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation).\n", + "Note that we need to make sure that all the edges of the polygons are also edges of the triangulation. We do this by using a [constrained (PSLG) Delaunay triangulation](https://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation).\n", "\n", "The result is a set of triangles covering the nested Grids, which we can use to determine in which Grid a particle is located at any given time. It is important that the list of polygons is ordered from smallest to largest Grid, so that triangles in overlapping areas are assigned to the correct Grid." ] @@ -242,7 +242,10 @@ " face_poly[ti] = ip\n", " break\n", "\n", - " return pts, tris, face_poly" + " return pts, tris, face_poly\n", + "\n", + "\n", + "points, face_tris, face_poly = constrained_triangulate_keep_edges(grid_polygons)" ] }, { @@ -250,7 +253,7 @@ "id": "15", "metadata": {}, "source": [ - "We can then run the triangulation and plot the resulting triangles to verify that they correctly cover the nested Grids." + "We can then plot the resulting triangles to verify that they correctly cover the nested Grids." ] }, { @@ -260,12 +263,10 @@ "metadata": {}, "outputs": [], "source": [ - "points, face_tris, face_poly = constrained_triangulate_keep_edges(grid_polygons)\n", - "\n", "fig, ax = plt.subplots()\n", "for i in range(n_grids)[::-1]:\n", " tris = face_tris[face_poly == i]\n", - " ax.triplot(points[:, 0], points[:, 1], tris, label=f\"Nest {i}\", color=f\"C{i}\")\n", + " ax.triplot(points[:, 0], points[:, 1], tris, label=f\"Grid {i}\", color=f\"C{i}\")\n", "ax.scatter(points[:, 0], points[:, 1], s=10, c=\"k\")\n", "ax.set_aspect(\"equal\")\n", "ax.legend()\n", @@ -281,7 +282,7 @@ "id": "17", "metadata": {}, "source": [ - "Then, we convert the triangulation into a Parcels FieldSet using `Parcels.UxGrid()`." + "Then, we convert the triangulation into an (unstructured) `parcels.FieldSet`." ] }, { @@ -309,7 +310,7 @@ " ),\n", " face_poly[np.newaxis, np.newaxis, :],\n", " {\n", - " \"long_name\": \"Nested grid ID\",\n", + " \"long_name\": \"Grid ID\",\n", " \"location\": \"face\",\n", " \"mesh\": \"delaunay\",\n", " },\n", @@ -326,14 +327,14 @@ "\n", "uxda = ux.UxDataArray(ds_tri[\"face_polygon\"], uxgrid=ux.Grid(ds_tri))\n", "\n", - "NestID = parcels.Field(\n", - " \"NestID\",\n", + "GridID = parcels.Field(\n", + " \"GridID\",\n", " uxda,\n", " # TODO note that here we need to use mesh=\"flat\" otherwise the hashing doesn't work. See https://github.com/Parcels-code/Parcels/pull/2439#discussion_r2627664010\n", " parcels.UxGrid(uxda.uxgrid, z=uxda[\"nz\"], mesh=\"flat\"),\n", " interp_method=parcels.interpolators.UxPiecewiseConstantFace,\n", ")\n", - "fieldset = parcels.FieldSet([NestID])" + "fieldset = parcels.FieldSet([GridID])" ] }, { @@ -341,7 +342,7 @@ "id": "19", "metadata": {}, "source": [ - "We can confirm that the FieldSet has been created correctly by running a Parcels simulation where particles sample the `NestID` field, which indicates in which nest each particle is located at any given time." + "We can confirm that the FieldSet has been created correctly by running a Parcels simulation where particles sample the `GridID` field, which indicates in which Grid each particle is located at any given time." ] }, { @@ -356,17 +357,24 @@ " np.linspace(-18, 38, 20),\n", ")\n", "\n", - "NestParticle = parcels.Particle.add_variable(parcels.Variable(\"nestID\", dtype=np.int32))\n", + "NestedGridParticle = parcels.Particle.add_variable(\n", + " parcels.Variable(\"gridID\", dtype=np.int32)\n", + ")\n", "pset = parcels.ParticleSet(\n", - " fieldset, pclass=NestParticle, lon=X.flatten(), lat=Y.flatten()\n", + " fieldset, pclass=NestedGridParticle, lon=X.flatten(), lat=Y.flatten()\n", ")\n", "\n", "\n", - "def SampleNestID(particles, fieldset):\n", - " particles.nestID = fieldset.NestID[particles]\n", + "def SampleGridID(particles, fieldset):\n", + " particles.gridID = fieldset.GridID[particles]\n", "\n", "\n", - "pset.execute(SampleNestID, runtime=np.timedelta64(1, \"s\"), dt=np.timedelta64(1, \"s\"))" + "pset.execute(\n", + " SampleGridID,\n", + " runtime=np.timedelta64(1, \"s\"),\n", + " dt=np.timedelta64(1, \"s\"),\n", + " verbose_progress=False,\n", + ")" ] }, { @@ -374,7 +382,7 @@ "id": "21", "metadata": {}, "source": [ - "Indeed, the visualisation below shows that particles correctly identify the nest they are in based on their location." + "Indeed, the visualisation below shows that particles correctly identify the grid they are in based on their location." ] }, { @@ -385,26 +393,21 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n", + "plot_args = {\"cmap\": cmap, \"edgecolors\": \"k\", \"linewidth\": 0.5}\n", "\n", "triang = mtri.Triangulation(\n", " uxda.uxgrid.node_lon.values,\n", " uxda.uxgrid.node_lat.values,\n", " triangles=uxda.uxgrid.face_node_connectivity.values,\n", ")\n", + "facecolors = np.squeeze(uxda[0, :].values)\n", "\n", - "plot_args = {\n", - " \"cmap\": cmap,\n", - " \"edgecolors\": \"k\",\n", - " \"linewidth\": 0.5,\n", - " \"vmin\": 0,\n", - " \"vmax\": 2.0,\n", - "}\n", - "ax.tripcolor(\n", - " triang, facecolors=np.squeeze(uxda[0, :].values), shading=\"flat\", **plot_args\n", - ")\n", - "ax.scatter(pset.lon, pset.lat, c=pset.nestID, **plot_args)\n", + "ax.tripcolor(triang, facecolors=facecolors, shading=\"flat\", **plot_args)\n", + "ax.scatter(pset.lon, pset.lat, c=pset.gridID, **plot_args)\n", "ax.set_aspect(\"equal\")\n", - "ax.set_title(\"Nesting visualisation (triangulation and interpolated particle values)\")\n", + "ax.set_title(\n", + " \"Nested Grids visualisation (triangulation and interpolated particle values)\"\n", + ")\n", "plt.tight_layout()\n", "plt.show()" ] @@ -414,11 +417,11 @@ "id": "23", "metadata": {}, "source": [ - "## Advecting particles with nest transitions\n", + "## Advecting particles with nested Grid transitions\n", "\n", - "We can now set up a particle advection simulation using the nested grids. We first combine all the Fields into a single FieldSet. \n", + "We can now set up a particle advection simulation using the nested Grids. We first combine all the Fields into a single FieldSet. \n", "\n", - "We rename the individual Fields by appending the nest index to their names, so that we can easily identify which Field belongs to which nest. We also add the `NestID` Field to the FieldSet (note that Parcels v4 supports combining structured and unstructured Fields into one FieldSet, which is very convenient for this usecase)." + "We rename the individual Fields by appending the Grid index to their names, so that we can easily identify which Field belongs to which Grid. We also add the `GridID` Field to the FieldSet (note that Parcels v4 supports combining structured and unstructured Fields into one FieldSet, which is very convenient for this usecase)." ] }, { @@ -428,14 +431,13 @@ "metadata": {}, "outputs": [], "source": [ - "fields = [NestID]\n", + "fields = [GridID]\n", "for i, ds in enumerate(ds_in):\n", " # TODO : use FieldSet.from_sgrid_convetion here once #2437 is merged\n", " grid = parcels.XGrid.from_dataset(ds, mesh=\"spherical\")\n", " U = parcels.Field(\"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.XLinear)\n", " V = parcels.Field(\"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.XLinear)\n", " UV = parcels.VectorField(\"UV\", U, V)\n", - "\n", " fset = parcels.FieldSet([U, V, UV])\n", "\n", " for fld in fset.fields.values():\n", @@ -449,7 +451,7 @@ "id": "25", "metadata": {}, "source": [ - "We then define a custom Advection kernel that advects particles using the appropriate velocity Field based on the `NestID` at the particle's location. Note that for simplicity, we use Eulerian advection here, but in a real-world application you would typically use a higher-order scheme." + "We then define a custom Advection kernel that advects particles using the appropriate velocity Field based on the `GridID` at the particle's location. Note that for simplicity, we use Eulerian advection here, but in a real-world application you would typically use a higher-order scheme." ] }, { @@ -459,8 +461,8 @@ "metadata": {}, "outputs": [], "source": [ - "def AdvectEE_Nests(particles, fieldset):\n", - " particles.nestID = fieldset.NestID[particles]\n", + "def AdvectEE_NestedGrids(particles, fieldset):\n", + " particles.gridID = fieldset.GridID[particles]\n", "\n", " # TODO because of KernelParticle bug (GH #2143), we need to copy lon/lat/time to local variables\n", " time = particles.time\n", @@ -470,10 +472,10 @@ " u = np.zeros_like(particles.lon)\n", " v = np.zeros_like(particles.lat)\n", "\n", - " unique_ids = np.unique(particles.nestID)\n", - " for nid in unique_ids:\n", - " mask = particles.nestID == nid\n", - " UVField = getattr(fieldset, f\"UV{nid}\")\n", + " unique_ids = np.unique(particles.gridID)\n", + " for gid in unique_ids:\n", + " mask = particles.gridID == gid\n", + " UVField = getattr(fieldset, f\"UV{gid}\")\n", " (u[mask], v[mask]) = UVField[time[mask], z[mask], lat[mask], lon[mask]]\n", "\n", " particles.dlon += u * particles.dt\n", @@ -487,18 +489,15 @@ " )\n", "\n", "\n", - "def DeleteErrorParticles(particles, fieldset):\n", - " any_error = particles.state >= 50 # This captures all Errors\n", - " particles[any_error].state = parcels.StatusCode.Delete\n", - "\n", - "\n", "lat = np.linspace(-17, 35, 10)\n", "lon = np.full(len(lat), -5)\n", "\n", - "pset = parcels.ParticleSet(fieldset, pclass=NestParticle, lon=lon, lat=lat)\n", - "ofile = parcels.ParticleFile(\"nest_particles.zarr\", outputdt=np.timedelta64(1, \"D\"))\n", + "pset = parcels.ParticleSet(fieldset, pclass=NestedGridParticle, lon=lon, lat=lat)\n", + "ofile = parcels.ParticleFile(\n", + " \"nestedgrid_particles.zarr\", outputdt=np.timedelta64(1, \"D\")\n", + ")\n", "pset.execute(\n", - " [AdvectEE_Nests, DeleteErrorParticles],\n", + " AdvectEE_NestedGrids,\n", " runtime=np.timedelta64(60, \"D\"),\n", " dt=np.timedelta64(1, \"D\"),\n", " output_file=ofile,\n", @@ -513,12 +512,12 @@ "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n", + "fig, ax = plt.subplots(1, 1, figsize=(10, 4))\n", "\n", - "ds_out = xr.open_zarr(\"nest_particles.zarr\")\n", + "ds_out = xr.open_zarr(\"nestedgrid_particles.zarr\")\n", "\n", "plt.plot(ds_out.lon.T, ds_out.lat.T, \"k\", linewidth=0.5)\n", - "sc = ax.scatter(ds_out.lon, ds_out.lat, c=ds_out.nestID, s=4, cmap=cmap, vmin=0, vmax=2)\n", + "sc = ax.scatter(ds_out.lon, ds_out.lat, c=ds_out.gridID, s=4, cmap=cmap, vmin=0, vmax=2)\n", "xl, yl = ax.get_xlim(), ax.get_ylim()\n", "\n", "for i in range(n_grids - 1):\n", @@ -534,8 +533,8 @@ "ax.set_aspect(\"equal\")\n", "\n", "cbar = plt.colorbar(sc, ticks=[0, 1, 2], ax=ax)\n", - "cbar.set_label(\"Nest ID\")\n", - "ax.set_title(\"Particle advection through nests\")\n", + "cbar.set_label(\"Grid ID\")\n", + "ax.set_title(\"Particle advection through nested Grids\")\n", "plt.tight_layout\n", "plt.show()" ] @@ -545,9 +544,9 @@ "id": "28", "metadata": {}, "source": [ - "Indeed, we can see that the particles follow the cosine oscillation pattern in the coarser grid, and move purely zonally in the finer grid.\n", + "Indeed, we can see that the particles follow the cosine oscillation pattern in the coarsest grid, move northeast in the finer grid, and move purely zonally in the finest grid.\n", "\n", - "Note that in the plot above the particles at higher latitudes move farther in longitude-space (at a constant zonal velocity) because of the curvature of the earth, so that is expected." + "Note that in the plot above the particles at higher latitudes move farther eastward because of the curvature of the earth, so that is expected." ] } ], From 054d7357ed139e128d030ee2884abb42b5bf5e7d Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 18 Dec 2025 08:19:47 +0100 Subject: [PATCH 09/12] Making one of the grids a pentagonal shape --- .../examples/tutorial_nestedgrids.ipynb | 23 ++++++++----------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index 40eb8315a..8c54d7ee7 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -7,12 +7,9 @@ "source": [ "# 🖥️ Nested Grids\n", "\n", - "In some applications, you may have access to fields on different grids that each cover only part of the region of interest. Then, you would like to combine them all together. You may also have a grid covering the entire region and another one only covering part of it, but with a higher resolution. The set of those grids form what we call nested grids.\n", + "In some applications, you may have access to fields on different grids that each cover only part of the region of interest. Then, you would like to combine them all together. You may for example have one grid covering the entire region and one or more others that cover only part of the region at a finer resolution. The set of those grids form what we call **nested grids**.\n", "\n", - "In Parcels v4, we can use the new `uxarray` integration to determine in which grid a particle is located. We will demonstrate how to set up a simulation with multiple nested grids, and how to handle particle transitions between these grids.\n", - "\n", - "\n", - "This tutorial shows how to use these Nested Fields with a very idealised example." + "In Parcels v4, we can use the new `uxarray` integration to determine in which grid a particle is located. We will demonstrate how to set up a simulation with multiple nested grids, and how to handle particle transitions between these grids, with a very idealised example." ] }, { @@ -69,7 +66,7 @@ "id": "4", "metadata": {}, "source": [ - "Now we create a zonal and meridional velocity field defined on a small grid. Both the zonal (1m/s) and meridional (0 m/s) velocity are uniform." + "Now we create a set of zonal and meridional velocities defined on a small grid. Both the zonal (1m/s) and meridional (0 m/s) velocity are uniform." ] }, { @@ -79,7 +76,7 @@ "metadata": {}, "outputs": [], "source": [ - "polygon = np.array([(10, 15), (25, 10), (25, 25), (10, 35)])\n", + "polygon = np.array([(10, 15), (25, 10), (25, 25), (17, 35), (10, 32)])\n", "grid_polygons.append(polygon)\n", "\n", "ds_in.append(setup_nested_ds(polygon))\n", @@ -92,7 +89,7 @@ "id": "6", "metadata": {}, "source": [ - "Then, we define another set of zonal and meridional velocities on a slightly larger grid. In this case, the meridional velocity is also 1 m/s." + "Then, we define another set of zonal and meridional velocities on a slightly larger grid. In this case, both the zonal and meridional velocity are 1 m/s." ] }, { @@ -115,7 +112,7 @@ "id": "8", "metadata": {}, "source": [ - "Finally, we define another set of velocities on an even larger grid. The zonal velocity is the same as for the smaller grids, but the meridional velocity is now a cosine as a function of longitude." + "Finally, we define another set of velocities on an even larger grid. The zonal velocity is again 1 m/s, but the meridional velocity is now a cosine as a function of longitude." ] }, { @@ -152,7 +149,7 @@ "n_grids = len(ds_in)\n", "cmap = ListedColormap([f\"C{i}\" for i in range(n_grids)])\n", "\n", - "fig, ax = plt.subplots()\n", + "fig, ax = plt.subplots(figsize=(10, 5))\n", "\n", "for i in range(n_grids):\n", " ds = ds_in[i].isel(time=0, depth=0)\n", @@ -174,7 +171,7 @@ "id": "12", "metadata": {}, "source": [ - "Note, as seen in the plot above, that the dataset domains in this case are rectangular, but the polygons that will later define the nested Grid boundaries don't have to be. So we can even use this method to subset parts of a Grid." + "Note, as seen in the plot above, that the dataset domains in this case are rectangular, but the polygons that will later define the nested Grid boundaries don't have to be. They can have any shape! This means that we can also use this method to subset parts of a Grid, or to define subregions of a domain." ] }, { @@ -263,7 +260,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots()\n", + "fig, ax = plt.subplots(figsize=(10, 5))\n", "for i in range(n_grids)[::-1]:\n", " tris = face_tris[face_poly == i]\n", " ax.triplot(points[:, 0], points[:, 1], tris, label=f\"Grid {i}\", color=f\"C{i}\")\n", @@ -512,7 +509,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots(1, 1, figsize=(10, 4))\n", + "fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n", "\n", "ds_out = xr.open_zarr(\"nestedgrid_particles.zarr\")\n", "\n", From a2706ce8e2b23de8de925a73f5cfa1c7e6fa71f9 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 18 Dec 2025 08:25:13 +0100 Subject: [PATCH 10/12] Using fieldset.GridID for triangles plot --- docs/user_guide/examples/tutorial_nestedgrids.ipynb | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index 8c54d7ee7..7cb04b333 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -76,7 +76,7 @@ "metadata": {}, "outputs": [], "source": [ - "polygon = np.array([(10, 15), (25, 10), (25, 25), (17, 35), (10, 32)])\n", + "polygon = np.array([(10, 15), (25, 10), (25, 25), (17, 36), (10, 32)])\n", "grid_polygons.append(polygon)\n", "\n", "ds_in.append(setup_nested_ds(polygon))\n", @@ -393,9 +393,9 @@ "plot_args = {\"cmap\": cmap, \"edgecolors\": \"k\", \"linewidth\": 0.5}\n", "\n", "triang = mtri.Triangulation(\n", - " uxda.uxgrid.node_lon.values,\n", - " uxda.uxgrid.node_lat.values,\n", - " triangles=uxda.uxgrid.face_node_connectivity.values,\n", + " fieldset.GridID.grid.uxgrid.node_lon.values,\n", + " fieldset.GridID.grid.uxgrid.node_lat.values,\n", + " triangles=fieldset.GridID.grid.uxgrid.face_node_connectivity.values,\n", ")\n", "facecolors = np.squeeze(uxda[0, :].values)\n", "\n", From 597c70aff57d0c225d958d940c653116001cedd8 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 18 Dec 2025 14:10:34 +0100 Subject: [PATCH 11/12] Removing NestedField tests and updating migration guide --- docs/user_guide/v4-migration.md | 1 + tests-v3/test_fieldset_sampling.py | 77 ------------------------------ 2 files changed, 1 insertion(+), 77 deletions(-) diff --git a/docs/user_guide/v4-migration.md b/docs/user_guide/v4-migration.md index a5685e58e..44633515b 100644 --- a/docs/user_guide/v4-migration.md +++ b/docs/user_guide/v4-migration.md @@ -39,6 +39,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat - `Field.eval()` returns an array of floats instead of a single float (related to the vectorization) - `Field.eval()` does not throw OutOfBounds or other errors +- The `NestedField` class has been removed. See the Nested Grids how-to guide for how to set up Nested Grids in v4. ## GridSet diff --git a/tests-v3/test_fieldset_sampling.py b/tests-v3/test_fieldset_sampling.py index 11dec1988..1a3c68623 100644 --- a/tests-v3/test_fieldset_sampling.py +++ b/tests-v3/test_fieldset_sampling.py @@ -14,7 +14,6 @@ Geographic, Particle, ParticleSet, - StatusCode, Variable, ) from tests.utils import create_fieldset_global @@ -796,82 +795,6 @@ def test_multiple_grid_addlater_error(): assert fail -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="Implementation of NestedFields is being reconsidered in v4.") -def test_nestedfields(): - from parcels import NestedField - - xdim = 10 - ydim = 20 - - U1 = Field( - "U1", - 0.1 * np.ones((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - V1 = Field( - "V1", - 0.2 * np.ones((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - U2 = Field( - "U2", - 0.3 * np.ones((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 2.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 2.0, ydim, dtype=np.float32), - ) - V2 = Field( - "V2", - 0.4 * np.ones((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 2.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 2.0, ydim, dtype=np.float32), - ) - U = NestedField("U", [U1, U2]) - V = NestedField("V", [V1, V2]) - fieldset = FieldSet(U, V) - - P1 = Field( - "P1", - 0.1 * np.ones((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - P2 = Field( - "P2", - 0.2 * np.ones((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 2.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 2.0, ydim, dtype=np.float32), - ) - P = NestedField("P", [P1, P2]) - fieldset.add_field(P) - - def Recover(particle, fieldset, time): # pragma: no cover - if particle.state == StatusCode.ErrorOutOfBounds: - particle.dlon = 0 - particle.dlat = 0 - particle.ddepth = 0 - particle.lon = 0 - particle.lat = 0 - particle.p = 999 - particle.state = StatusCode.Evaluate - - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0], lat=[0.3]) - pset.execute(AdvectionRK4 + pset.Kernel(SampleP), runtime=2, dt=1) - assert np.isclose(pset.lat[0], 0.5) - assert np.isclose(pset.p[0], 0.1) - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0], lat=[1.1]) - pset.execute(AdvectionRK4 + pset.Kernel(SampleP), runtime=2, dt=1) - assert np.isclose(pset.lat[0], 1.5) - assert np.isclose(pset.p[0], 0.2) - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0], lat=[2.3]) - pset.execute(pset.Kernel(AdvectionRK4) + SampleP + Recover, runtime=1, dt=1) - assert np.isclose(pset.lat[0], 0) - assert np.isclose(pset.p[0], 999) - assert np.allclose(fieldset.UV[0][0, 0, 0, 0], [0.1, 0.2]) - - def test_fieldset_sampling_updating_order(tmp_zarrfile): def calc_p(t, y, x): return 10 * t + x + 0.2 * y From 3d1be3e24fcf877d3127c2c7cf2354e150e7b58b Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 19 Dec 2025 12:53:07 +0100 Subject: [PATCH 12/12] Implementing review suggestions Review suggestions by @reint-fischer --- .../examples/tutorial_nestedgrids.ipynb | 66 +++++++++++-------- 1 file changed, 38 insertions(+), 28 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nestedgrids.ipynb b/docs/user_guide/examples/tutorial_nestedgrids.ipynb index 7cb04b333..7cc465100 100644 --- a/docs/user_guide/examples/tutorial_nestedgrids.ipynb +++ b/docs/user_guide/examples/tutorial_nestedgrids.ipynb @@ -42,10 +42,20 @@ "First we define a helper function to quickly set up a nested dataset" ] }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "```{note}\n", + "If your nested data come from hydrodynamic models, you only have to provide the `grid_polygons` list with polygon coordinates for each grid, but don't need to create the `ds_in` like below.\n", + "```" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "3", + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -63,7 +73,7 @@ }, { "cell_type": "markdown", - "id": "4", + "id": "5", "metadata": {}, "source": [ "Now we create a set of zonal and meridional velocities defined on a small grid. Both the zonal (1m/s) and meridional (0 m/s) velocity are uniform." @@ -72,7 +82,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -86,7 +96,7 @@ }, { "cell_type": "markdown", - "id": "6", + "id": "7", "metadata": {}, "source": [ "Then, we define another set of zonal and meridional velocities on a slightly larger grid. In this case, both the zonal and meridional velocity are 1 m/s." @@ -95,7 +105,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -109,7 +119,7 @@ }, { "cell_type": "markdown", - "id": "8", + "id": "9", "metadata": {}, "source": [ "Finally, we define another set of velocities on an even larger grid. The zonal velocity is again 1 m/s, but the meridional velocity is now a cosine as a function of longitude." @@ -118,7 +128,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -133,7 +143,7 @@ }, { "cell_type": "markdown", - "id": "10", + "id": "11", "metadata": {}, "source": [ "We plot the velocity fields on top of each other, indicating the grid boundaries in red." @@ -142,7 +152,7 @@ { "cell_type": "code", "execution_count": null, - "id": "11", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -159,7 +169,7 @@ " np.append(poly[:, 0], poly[0, 0]),\n", " np.append(poly[:, 1], poly[0, 1]),\n", " \"-r\",\n", - " lw=1,\n", + " lw=2,\n", " )\n", "\n", "plt.tight_layout()\n", @@ -168,7 +178,7 @@ }, { "cell_type": "markdown", - "id": "12", + "id": "13", "metadata": {}, "source": [ "Note, as seen in the plot above, that the dataset domains in this case are rectangular, but the polygons that will later define the nested Grid boundaries don't have to be. They can have any shape! This means that we can also use this method to subset parts of a Grid, or to define subregions of a domain." @@ -176,12 +186,12 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "14", "metadata": {}, "source": [ "## Creating a Delaunay triangulation of the nested Grids\n", "\n", - "Now comes the important part: we need to create a Delaunay triangulation of the nested Grids, so that we can efficiently determine in which Grid a particle is located at any given time. We use the `triangle` package to perform the triangulation, and `shapely` to handle the geometric operations. \n", + "Now comes the important part: we need to create a Delaunay triangulation of the nested Grids, so that we can efficiently determine in which Grid a particle is located at any given time. We use the [`triangle` package](https://rufat.be/triangle/delaunay.html) to perform the triangulation, and `shapely` to handle the geometric operations. \n", "\n", "Note that we need to make sure that all the edges of the polygons are also edges of the triangulation. We do this by using a [constrained (PSLG) Delaunay triangulation](https://en.wikipedia.org/wiki/Constrained_Delaunay_triangulation).\n", "\n", @@ -191,7 +201,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "15", "metadata": {}, "outputs": [], "source": [ @@ -247,7 +257,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "16", "metadata": {}, "source": [ "We can then plot the resulting triangles to verify that they correctly cover the nested Grids." @@ -256,7 +266,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "17", "metadata": {}, "outputs": [], "source": [ @@ -276,7 +286,7 @@ }, { "cell_type": "markdown", - "id": "17", + "id": "18", "metadata": {}, "source": [ "Then, we convert the triangulation into an (unstructured) `parcels.FieldSet`." @@ -285,7 +295,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "19", "metadata": {}, "outputs": [], "source": [ @@ -336,7 +346,7 @@ }, { "cell_type": "markdown", - "id": "19", + "id": "20", "metadata": {}, "source": [ "We can confirm that the FieldSet has been created correctly by running a Parcels simulation where particles sample the `GridID` field, which indicates in which Grid each particle is located at any given time." @@ -345,7 +355,7 @@ { "cell_type": "code", "execution_count": null, - "id": "20", + "id": "21", "metadata": {}, "outputs": [], "source": [ @@ -376,7 +386,7 @@ }, { "cell_type": "markdown", - "id": "21", + "id": "22", "metadata": {}, "source": [ "Indeed, the visualisation below shows that particles correctly identify the grid they are in based on their location." @@ -385,7 +395,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22", + "id": "23", "metadata": {}, "outputs": [], "source": [ @@ -411,7 +421,7 @@ }, { "cell_type": "markdown", - "id": "23", + "id": "24", "metadata": {}, "source": [ "## Advecting particles with nested Grid transitions\n", @@ -424,7 +434,7 @@ { "cell_type": "code", "execution_count": null, - "id": "24", + "id": "25", "metadata": {}, "outputs": [], "source": [ @@ -445,7 +455,7 @@ }, { "cell_type": "markdown", - "id": "25", + "id": "26", "metadata": {}, "source": [ "We then define a custom Advection kernel that advects particles using the appropriate velocity Field based on the `GridID` at the particle's location. Note that for simplicity, we use Eulerian advection here, but in a real-world application you would typically use a higher-order scheme." @@ -454,7 +464,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26", + "id": "27", "metadata": {}, "outputs": [], "source": [ @@ -505,7 +515,7 @@ { "cell_type": "code", "execution_count": null, - "id": "27", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -538,7 +548,7 @@ }, { "cell_type": "markdown", - "id": "28", + "id": "29", "metadata": {}, "source": [ "Indeed, we can see that the particles follow the cosine oscillation pattern in the coarsest grid, move northeast in the finer grid, and move purely zonally in the finest grid.\n",