Cleaning up text and code on nestedgrid tutorial

Author: erikvansebille

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,15 +242,18 @@
     "                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)"
    ]
   },
   {
    "cell_type": "markdown",
    "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,22 +327,22 @@
     "\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])"
    ]
   },
   {
    "cell_type": "markdown",
    "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,25 +357,32 @@
     "    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",
+    ")"
    ]
   },
   {
    "cell_type": "markdown",
    "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."
    ]
   }
  ],