Grid explanations user guide (as ipynb)

docs/user_guide/examples/explanation_grids.ipynb

@@ -0,0 +1,345 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# 📖 Grids\n",
+    "\n",
+    "In Lagrangian ocean analysis, virtual particle tracking requires the accurate interpolation of physical properties (flow velocities and tracer properties) to particle locations. The underlying data that forces the particle movement will likely be defined on a discretised grid. \n",
+    "\n",
+    "Parcels can natively handle two styles of grids; structured and unstructured, where Parcels `Field` objects exist on a (structured) `parcels.XGrid` and conform to [SGRID](https://sgrid.github.io/sgrid/) conventions, or on a (unstructured) `parcels.Uxgrid` and conform to [UGRID](https://ugrid-conventions.github.io/ugrid-conventions/) conventions. Here we describe these grids on a conceptual level.\n",
+    "\n",
+    "```{note}\n",
+    "Parcels comes ready with a range of `convert` functions which will convert your `xarray` dataset into a SGRID or UGRID compliant dataset. See [this list](../index.md#converting-model-data-to-fieldsets) for all the models we have `convert` functions for. If your model is not on that list, you can create your own `convert` function, or reach out on our discussion board for help.\n",
+    "```\n",
+    "\n",
+    "Under the hood, every `Field` in a `FieldSet` has a `grid` attribute. This `grid` stores the spatial and temporal information of the Field coordinates. The number of Grids in a FieldSet is thus always smaller or equal to the number of Field objects; and this is what the \"grid number\" column in `FieldSet.describe()` refers to.\n",
+    "\n",
+    "## Structured grids\n",
+    "\n",
+    "A structured grid is composed of quadrilateral elements, that are indexed using 3D indices, like $(i,j,k)$. In `xarray` terminology, these indices correspond to [`dimensions`](https://docs.xarray.dev/en/latest/user-guide/terminology.html#term-Dimension). For example, in NEMO datasets, $(x,y)$ typically define these dimensions horizontally, where the physical coordinates are $(glamf, gphif)$. \n",
+    "\n",
+    "A major benefit of structured grids is that grid cell neighbours are easily found by decrementing or incrementing these dimensions. However, in Parcels we either perform a binary search in the case of 1-dimensional coordinates, or use a [hash table](../../development/unstructured_grid_search.md) in the case of 2/3-dimensional coordinates.\n",
+    "\n",
+    "There are two styles of structured grids, rectilinear and curvilinear, as shown in Figure 1below.\n",
+    "\n",
+    "1. Rectilinear grids are typically aligned with the coordinate axes (that is, there is a one-to-one mapping between the dimensions and the physical coordinate space), making them simple and computationally efficient to query. However, it is difficult to resolve complex coastlines without high resolution, and they suffer from singularities in the flow fields at the poles where the grid lines converge.\n",
+    "\n",
+    "2. Curvilinear grids allow for curved grid lines, which allow for better representation of coastlines. Poles are often on land to avoid singularities in the ocean. As the spatial positions of their vertices (their [\"coordinates\" in xarray-parlance](https://docs.xarray.dev/en/latest/user-guide/terminology.html#term-Coordinate)) are stored in separate 2D arrays, moving to an eastward neighbour is not as simple as incrementing the $i$-th dimension by 1.\n",
+    "\n",
+    "<figure>\n",
+    "  <img src=\"../../_static/grid_discretisation.png\" alt=\"Grid discretisations handled by Parcels\" />\n",
+    "  <figcaption>Figure 1 - Grid discretisations handled by Parcels. In the horizontal plane; (a) rectilinear, (b) curvilinear. In the vertical plane; (c) z-levels, (d) sigma-levels. Adapted from the <a href=\"https://doi.org/10.5194/gmd-12-3571-2019\">Parcels v2.0 paper</a></figcaption>\n",
+    "</figure>\n",
+    "\n",
+    "## Unstructured grids\n",
+    "\n",
+    "TODO: For Joe.\n",
+    "\n",
+    "## How your data may be defined\n",
+    "\n",
+    "Regardless of your grid type, how your data is defined on your grid is equally important. Here, we will describe how you can interpret your data at a very general level, and the concept applies for any $n$-gon ($n$-sided polygon with $n \\ge 3$).\n",
+    "\n",
+    "### Nodes, edges, and faces\n",
+    "\n",
+    "Let's assume we have a simple 2D rectilinear grid. The grid is composed of a number of nodes (or vertices), connected by edges, which together construct grid cells and grid cell faces. In Figure 2, we draw a simple grid cell. Here, your data may be defined on the corners of the grid cell (the nodes/vertices), at the centre of the cell face, or across a cell edge."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# Ignore RuntimeWarnings that arise in plotting\n",
+    "import warnings\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "from scipy.interpolate import RegularGridInterpolator, interp1d\n",
+    "\n",
+    "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n",
+    "\n",
+    "plt.figure(figsize=(4, 4))\n",
+    "plt.axvline(x=0, color=\"black\")\n",
+    "plt.axvline(x=1, color=\"black\")\n",
+    "plt.axhline(y=0, color=\"black\")\n",
+    "plt.axhline(y=1, color=\"black\")\n",
+    "plt.scatter([0, 0, 1, 1], [0, 1, 0, 1], color=\"blue\", s=50, label=\"Nodes/Vertices\")\n",
+    "plt.scatter([0.5], [0.5], color=\"red\", s=50, label=\"Centre of grid cell\")\n",
+    "plt.xlim(-0.25, 1.25)\n",
+    "plt.ylim(-0.25, 1.25)\n",
+    "plt.xticks([])\n",
+    "plt.yticks([])\n",
+    "plt.legend()\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2",
+   "metadata": {},
+   "source": [
+    "Figure 2 - A simple grid cell. Blue circles denote nodes (or vertices) of the grid cell. The red circle denotes the centre of the cell face.\n",
+    "\n",
+    "You will need to make several assumptions about your data. Your data may represent a point-wise \"sample\" of some field. For example, velocity data may be defined at the nodes of your grid, and a typical assumption to make is that you can bi-linearly interpolate these data points to your particle positions. In such an example, your velocity field may look like Figure 3. Bi-linear interpolation ensures continuity of the velocity at the cell boundaries, however, it does not ensure a smooth (differentiable) transition.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3",
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# Construct dummy data for velocity field at the nodes of a single grid cell\n",
+    "xcoords, ycoords = np.linspace(0, 1, 2), np.linspace(0, 1, 2)\n",
+    "\n",
+    "u_data = np.array([[0, -0.5], [1, -1]])\n",
+    "v_data = np.array([[0, -2], [0, -2]])\n",
+    "\n",
+    "interp_u = RegularGridInterpolator((xcoords, ycoords), u_data, method=\"linear\")\n",
+    "\n",
+    "interp_v = RegularGridInterpolator((xcoords, ycoords), v_data, method=\"linear\")\n",
+    "\n",
+    "interp_xx, interp_yy = np.linspace(0, 1, 51), np.linspace(0, 1, 51)\n",
+    "interp_xx, interp_yy = np.meshgrid(interp_xx, interp_yy, indexing=\"ij\")\n",
+    "\n",
+    "eval_points = np.stack([interp_xx.ravel(), interp_yy.ravel()], axis=-1)\n",
+    "\n",
+    "interp_u_data = interp_u(eval_points).reshape(51, 51)\n",
+    "interp_v_data = interp_v(eval_points).reshape(51, 51)\n",
+    "\n",
+    "interp_speed_data = np.sqrt(interp_u_data**2 + interp_v_data**2)\n",
+    "\n",
+    "# Plot the interpolated velocity vectors and speed data\n",
+    "mm = 5\n",
+    "plt.axhline(0, linewidth=0.5, color=\"k\")\n",
+    "plt.axvline(0, linewidth=0.5, color=\"k\")\n",
+    "plt.axhline(1, linewidth=0.5, color=\"k\")\n",
+    "plt.axvline(1, linewidth=0.5, color=\"k\")\n",
+    "cb = plt.pcolormesh(\n",
+    "    interp_xx, interp_yy, interp_speed_data, shading=\"gouraud\", cmap=plt.cm.viridis\n",
+    ")\n",
+    "plt.quiver(\n",
+    "    interp_xx[::mm, ::mm],\n",
+    "    interp_yy[::mm, ::mm],\n",
+    "    interp_u_data[::mm, ::mm],\n",
+    "    interp_v_data[::mm, ::mm],\n",
+    ")\n",
+    "plt.title(\"Bi-linear interpolation of velocity data at nodes\")\n",
+    "\n",
+    "plt.colorbar(cb, ax=plt.gca(), label=\"Speed [m/s]\")\n",
+    "plt.xlabel(\"X [km]\")\n",
+    "plt.ylabel(\"Y [km]\")\n",
+    "\n",
+    "plt.xlim([-0.2, 1.2])\n",
+    "plt.ylim([-0.2, 1.2])\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4",
+   "metadata": {},
+   "source": [
+    "**Figure 3 - Bi-linear interpolation of point-wise data at nodes.**\n",
+    "\n",
+    "Alternatively, your data may represent an \"average value\" across a cell face. For example, your temperature and salinity data may be defined at the cell centre, and represent an average value for the entire grid cell. A typical assumption to make is that you can nearest-neighbour interpolate these data points to your particle positions. In such a case, your temperature field may look like figure 4. A nearest-neighbour interpolation scheme ensures you have a piece-wise constant field, typically with sharp transitions at grid cell boundaries."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5",
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "xcoords, ycoords = np.linspace(0, 1, 2), np.linspace(0, 1, 2)\n",
+    "xx, yy = np.meshgrid(xcoords, ycoords, indexing=\"ij\")\n",
+    "data = 18 + np.random.rand(1) * np.array([[-1, 1], [0, 2]])\n",
+    "\n",
+    "mm = 5\n",
+    "# Paste the data at different locations to show the effect of nearest-neighbour interpolation\n",
+    "for i in [1, -1]:\n",
+    "    for j in [-1, 1]:\n",
+    "        plt.pcolormesh(\n",
+    "            xx + i * 0.5, yy + j * 0.5, data, shading=\"auto\", cmap=plt.cm.viridis\n",
+    "        )\n",
+    "        plt.scatter(xx + i * 0.5, yy + j * 0.5, c=\"k\", zorder=20)\n",
+    "cb = plt.pcolormesh(\n",
+    "    xx + 0.5, yy + 0.5, 18 + 0.5 * data, shading=\"auto\", cmap=plt.cm.viridis, zorder=10\n",
+    ")\n",
+    "\n",
+    "plt.title(\"Nearest-neighbour interpolation of\\ntemperature data at cell centres\")\n",
+    "\n",
+    "plt.colorbar(cb, ax=plt.gca(), label=\"Temperature [deg C]\")\n",
+    "plt.xlabel(\"X [km]\")\n",
+    "plt.ylabel(\"Y [km]\")\n",
+    "plt.axhline(-1, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "plt.axhline(0, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "plt.axhline(1, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "plt.axhline(2, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "plt.axvline(-1, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "plt.axvline(0, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "plt.axvline(1, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "plt.axvline(2, linewidth=0.5, color=\"k\", zorder=30)\n",
+    "\n",
+    "plt.xlim([-1, 2])\n",
+    "plt.ylim([-1, 2])\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "metadata": {},
+   "source": [
+    "**Figure 4 - Nearest-neighbour interpolation of \"grid cell averaged\" data at cell centres.**\n",
+    "\n",
+    "Your data may represent a value across a cell edge. For example, in (2D) Arakawa C-grid datasets, velocities are defined across an edge as they represent a \"flux\" across that cell edge. Additionally, these cell edges may not be aligned with the coordinate axes, and rather represent a velocity in the $i$ or $j$ direction.\n",
+    "\n",
+    "For structured grids, [Blanke and Raynaud](<https://doi.org/10.1175/1520-0485(1997)027%3C1038:KOTPEU%3E2.0.CO;2>) proposed in 1997 to perform a 1D linear interpolation of the $i$ velocity in the $i$ direction, and similarly a 1D linear interpolation of the $j$ velocity in the $j$ direction. These velocities must then be rotated into meridional and zonal velocities, which Parcels handles under the hood. In such a case, your velocity field may look like figure 5.\n",
+    "\n",
+    "This (uni)linear velocity interpolation is now often referred to as the Analytical interpolation scheme, and is what the [Ariane](https://ariane-code.cnrs.fr) and [TRACMASS](https://www.tracmass.org/index.html) Lagrangian codes also use. In Parcels, the time-stepping version of this interpolation is provided in the `CGrid_Velocity` Interpolator function.\n",
+    "\n",
+    "```{note}\n",
+    "If you have a structured curvilinear grid, and your velocity field is oriented in the $i/j$-directions, then it is almost certain that your velocities are computed from fluxes across the cell edge, and that you will require the `CGRID_Velocity` Interpolator.\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7",
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# Construct dummy data for velocity field at the cell edges of a single grid cell\n",
+    "u_data = np.array([[2, 1]])\n",
+    "v_data = np.array([[0, 1]])\n",
+    "\n",
+    "interp_u = interp1d([0, 1], u_data.flatten(), kind=\"linear\")\n",
+    "interp_v = interp1d([0, 1], v_data.flatten(), kind=\"linear\")\n",
+    "\n",
+    "interp_xx, interp_yy = np.linspace(0, 1, 51), np.linspace(0, 1, 51)\n",
+    "XX, YY = np.meshgrid(interp_xx, interp_yy, indexing=\"xy\")\n",
+    "\n",
+    "interp_u_data = interp_u(XX)\n",
+    "interp_v_data = interp_v(YY)\n",
+    "interp_speed_data = np.sqrt(interp_u_data**2 + interp_v_data**2)\n",
+    "\n",
+    "# Plot the interpolated velocity vectors and speed data\n",
+    "mm = 5\n",
+    "plt.axhline(0, linewidth=0.5, color=\"k\")\n",
+    "plt.axvline(0, linewidth=0.5, color=\"k\")\n",
+    "plt.axhline(1, linewidth=0.5, color=\"k\")\n",
+    "plt.axvline(1, linewidth=0.5, color=\"k\")\n",
+    "\n",
+    "cb = plt.pcolormesh(interp_xx, interp_yy, interp_speed_data, vmin=0.75, vmax=2)\n",
+    "plt.quiver(\n",
+    "    interp_xx[::mm],\n",
+    "    interp_yy[::mm],\n",
+    "    interp_u_data[::mm, ::mm],\n",
+    "    interp_v_data[::mm, ::mm],\n",
+    ")\n",
+    "plt.quiver(\n",
+    "    interp_xx[0],\n",
+    "    interp_yy[len(interp_yy) // 2],\n",
+    "    interp_u_data[len(interp_u_data) // 2, 0],\n",
+    "    0,\n",
+    "    color=\"r\",\n",
+    "    scale=30,\n",
+    ")\n",
+    "plt.quiver(\n",
+    "    interp_xx[-1],\n",
+    "    interp_yy[len(interp_yy) // 2],\n",
+    "    interp_u_data[len(interp_u_data) // 2, -1],\n",
+    "    0,\n",
+    "    color=\"r\",\n",
+    "    scale=30,\n",
+    ")\n",
+    "\n",
+    "plt.quiver(\n",
+    "    interp_xx[len(interp_yy) // 2],\n",
+    "    interp_yy[0],\n",
+    "    0,\n",
+    "    interp_v_data[0, len(interp_v_data) // 2],\n",
+    "    color=\"b\",\n",
+    ")\n",
+    "plt.quiver(\n",
+    "    interp_xx[len(interp_yy) // 2],\n",
+    "    interp_yy[-1],\n",
+    "    0,\n",
+    "    interp_v_data[-1, len(interp_v_data) // 2],\n",
+    "    color=\"b\",\n",
+    ")\n",
+    "\n",
+    "plt.title(\"C-grid interpolation of velocity data\\nacross cell edges\")\n",
+    "plt.colorbar(cb, ax=plt.gca(), label=\"Speed [m/s]\")\n",
+    "plt.xlim([-0.2, 1.2])\n",
+    "plt.ylim([-0.2, 1.2])\n",
+    "plt.xlabel(\"X [km]\")\n",
+    "plt.ylabel(\"Y [km]\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8",
+   "metadata": {},
+   "source": [
+    "**Figure 5 - C-grid 1D interpolation of velocity data defined across a cell edge.**\n",
+    "\n",
+    "\n",
+    "### Vertical coordinates\n",
+    "\n",
+    "Lastly, a short note on vertical coordinates. Parcels can handle two styles of vertical coordinates; z-levels which define fixed depth levels in physical space, and sigma-levels which define varying depth levels in physical space as a function of the water column depth. As sigma-levels are effectively \"terrain-following\", the grid cell faces may no longer be orthogonal to the domain surface. Figure 1 visualises these differences, and parcels handles this all under the hood. See the CROCO Tutorial for more details on the sigma-levels implementation.\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "parcels_copernicus",
+   "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.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/user_guide/index.md

@@ -30,18 +30,26 @@ getting_started/tutorial_output.ipynb
 getting_started/explanation_concepts.md
 ```
 
-## Set up FieldSets
+## Converting model data to FieldSets
 
 ```{toctree}
-:caption: Set up FieldSets
-:name: setup-fieldsets
+:caption: Converting model data to FieldSets
+:name: converting-model-data
 :titlesonly:
-examples/explanation_grids.md
 examples/tutorial_nemo.ipynb
 examples/tutorial_croco_3D.ipynb
 examples/tutorial_mitgcm.ipynb
 examples/tutorial_fesom.ipynb
 examples/tutorial_schism.ipynb
+```
+
+## Work with FieldSets
+
+```{toctree}
+:caption: Work with FieldSets
+:name: work-with-fieldsets
+:titlesonly:
+examples/explanation_grids.ipynb
 examples/tutorial_velocityconversion.ipynb
 examples/tutorial_nestedgrids.ipynb
 examples/tutorial_manipulating_field_data.ipynb