diff --git a/docs/user_guide/examples/tutorial_schism.ipynb b/docs/user_guide/examples/tutorial_schism.ipynb new file mode 100644 index 000000000..432c98e7d --- /dev/null +++ b/docs/user_guide/examples/tutorial_schism.ipynb @@ -0,0 +1,435 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# 🖥️ SCHISM tutorial" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "Parcels v4 supports unstructured-grid model output via [uxarray](https://uxarray.readthedocs.io/).\n", + "This tutorial walks through advecting particles in real [SCHISM](http://ccrm.vims.edu/schismweb/)\n", + "output from a hydrodynamic hindcast of **Lake Ontario**. SCHISM writes its output in the\n", + "[scribe-IO](https://schism-dev.github.io/schism/master/getting-started/output.html) format, which\n", + "already follows the [UGRID](https://ugrid-conventions.github.io/ugrid-conventions/) conventions for\n", + "its horizontal mesh, so `uxarray` can read it directly.\n", + "\n", + "SCHISM output differs from the [FESOM tutorial](./tutorial_fesom.ipynb) in a few ways that this\n", + "tutorial highlights:\n", + "\n", + "1. The mesh coordinates are in a **projected coordinate system (metres)**, not longitude/latitude, so\n", + " we build the grid as a *flat* (Cartesian) mesh.\n", + "2. Velocities are stored at **mesh nodes** over a vertical column of layers, ordered **bottom → surface**.\n", + "3. SCHISM uses a **localized vertical grid (LSC2)**: the number of valid vertical levels varies from\n", + " node to node, and levels below the seabed are stored as `NaN`. We use this to flag particles that\n", + " drift below the local bathymetry and stop advecting them.\n", + "\n", + "If you have not done so already, work through the\n", + "[quickstart tutorial](../../getting_started/tutorial_quickstart.md) first to get familiar with\n", + "`ParticleSet`, `Kernel`, and `ParticleFile`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "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", + "\n", + "import parcels\n", + "import parcels.tutorial\n", + "from parcels._core.statuscodes import StatusCode" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## Get the SCHISM tutorial dataset\n", + "\n", + "We use three files from a SCHISM Lake Ontario hindcast (6 hourly snapshots), bundled in Parcels' tutorial\n", + "data registry:\n", + "\n", + "* `out2d`: the 2D output, slimmed to the **horizontal mesh** (node/face topology) and bathymetry.\n", + "* `horizontalVelX`, `horizontalVelY`: the 3D horizontal velocity components, defined at the mesh nodes\n", + " over 32 vertical layers.\n", + "\n", + "As in the [quickstart](../../getting_started/tutorial_quickstart.md), `parcels.tutorial.open_dataset`\n", + "downloads the files into a local cache on first use (subsequent calls return the cached copy) and opens\n", + "them as `xarray` datasets:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "grid_ds = parcels.tutorial.open_dataset(\"SCHISM_LakeOntario/out2d\")\n", + "u = parcels.tutorial.open_dataset(\"SCHISM_LakeOntario/horizontalVelX\")[\"horizontalVelX\"]\n", + "v = parcels.tutorial.open_dataset(\"SCHISM_LakeOntario/horizontalVelY\")[\"horizontalVelY\"]\n", + "grid_ds" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## Build the horizontal mesh\n", + "\n", + "SCHISM stores its mesh topology in `out2d_1.nc` following the UGRID conventions: node coordinates\n", + "(`SCHISM_hgrid_node_x/y`) and a face–node connectivity table (`SCHISM_hgrid_face_nodes`). Two\n", + "SCHISM-specific details need handling before we hand the mesh to Parcels:\n", + "\n", + "* **Triangular cells.** The connectivity table is stored with a width of 4 (so the same format can\n", + " describe quads), but this Lake Ontario mesh is entirely triangular; the 4th column is all fill. We\n", + " keep the first three columns and convert the 1-based indices to 0-based. Parcels' `UxGrid` requires\n", + " purely triangular cells.\n", + "* **Projected coordinates.** The coordinates are in metres (`standard_name = projection_x_coordinate`),\n", + " not degrees. `uxarray` currently assumes node coordinates are spherical and wraps longitudes into\n", + " [-180, 180] (see [uxarray #1524](https://github.com/UXARRAY/uxarray/issues/1524)), which would corrupt\n", + " the mesh. We undo that wrap by writing the raw metre coordinates back, and later build the `FieldSet`\n", + " with `mesh=\"flat\"` so Parcels treats the plane as Cartesian." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "node_x = grid_ds[\"SCHISM_hgrid_node_x\"].values.astype(\"float64\")\n", + "node_y = grid_ds[\"SCHISM_hgrid_node_y\"].values.astype(\"float64\")\n", + "face_nodes = (\n", + " grid_ds[\"SCHISM_hgrid_face_nodes\"].values[:, :3].astype(\"int64\") - 1\n", + ") # all-triangular, 0-based\n", + "\n", + "uxgrid = ux.Grid.from_topology(\n", + " node_lon=node_x, node_lat=node_y, face_node_connectivity=face_nodes, fill_value=-1\n", + ")\n", + "# undo uxarray's [-180, 180] longitude wrap of the projected metres (uxarray #1524)\n", + "uxgrid.node_lon.values[:] = node_x\n", + "uxgrid.node_lat.values[:] = node_y\n", + "\n", + "print(\n", + " f\"n_node={uxgrid.n_node}, n_face={uxgrid.n_face}, n_max_face_nodes={uxgrid.n_max_face_nodes}\"\n", + ")\n", + "print(f\"x range: {node_x.min():.0f} .. {node_x.max():.0f} m\")" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Assemble the velocity fields\n", + "\n", + "The two velocity files hold `horizontalVelX` and `horizontalVelY` with dimensions\n", + "`(time, node, layer)`. We rename them to `U`/`V` (so Parcels recognises the velocity components) and\n", + "to the Parcels UGRID dimension names (`n_node` for the lateral dimension).\n", + "\n", + "Two vertical details:\n", + "\n", + "* **Layer ordering.** SCHISM stores levels **bottom → surface**, while Parcels expects depth increasing\n", + " downward from the surface, so we reverse the layer axis.\n", + "* **Vertical coordinate.** SCHISM's LSC2 vertical grid (*Localized Sigma Coordinates with Shaved\n", + " cells*) varies with horizontal position: each node has its own layer depths, and even its own\n", + " *number* of levels (the true depths are written to a separate `zCoordinates` output, not used here).\n", + " **Parcels does not currently support a vertical grid that varies with lateral position.** `UxGrid`\n", + " takes a single 1D column of layer-interface depths that applies everywhere on the mesh. We therefore\n", + " supply a fictitious 1D vertical grid. This is adequate for the near-surface, horizontal transport\n", + " shown here (the lateral interpolation does not depend on the vertical grid); accurate full-depth 3D\n", + " transport on an LSC2 grid would require Parcels to support a laterally varying vertical coordinate.\n", + "\n", + "We also call `.load()` so the velocities sit in memory; otherwise every interpolation step re-reads\n", + "from disk and the simulation is extremely slow." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "nlev = u.sizes[\"nSCHISM_vgrid_layers\"]\n", + "\n", + "# placeholder interface depths (metres, positive down): surface (0 m) -> bed, refined near the surface\n", + "zf = 250.0 * (np.arange(nlev) / (nlev - 1)) ** 1.5\n", + "zc = 0.5 * (zf[1:] + zf[:-1])\n", + "\n", + "rename = {\"nSCHISM_vgrid_layers\": \"zf\", \"nSCHISM_hgrid_node\": \"n_node\"}\n", + "U = (\n", + " u.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename).load()\n", + ") # reverse to surface-first\n", + "V = v.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename).load()\n", + "\n", + "uxds = ux.UxDataset(\n", + " xr.Dataset(\n", + " {\n", + " \"U\": U.transpose(\"time\", \"zf\", \"n_node\"),\n", + " \"V\": V.transpose(\"time\", \"zf\", \"n_node\"),\n", + " },\n", + " coords={\"zf\": (\"zf\", zf), \"zc\": (\"zc\", zc)},\n", + " ),\n", + " uxgrid=uxgrid,\n", + ")\n", + "uxds" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Build the `FieldSet`\n", + "\n", + "With the mesh and UGRID-compliant dimensions in place, `parcels.FieldSet.from_ugrid_conventions` builds\n", + "the `FieldSet`. It detects `U` and `V`, attaches the `UxGrid`, and selects the `UxLinearNodeLinearZF`\n", + "interpolator (barycentric in the horizontal, linear in the vertical) because the velocities are\n", + "node-registered along the layer interfaces `zf`. We pass `mesh=\"flat\"` because the coordinates\n", + "are projected metres: velocities in m/s then advect positions in metres directly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "fieldset = parcels.FieldSet.from_ugrid_conventions(uxds, mesh=\"flat\")\n", + "\n", + "for name, field in fieldset.fields.items():\n", + " interp = getattr(field, \"interp_method\", None)\n", + " print(\n", + " f\"{name:>4s} -> {type(field).__name__:<11s} interp={interp.__name__ if interp else '-'}\"\n", + " )\n", + "print(\"time interval:\", fieldset.time_interval)" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Stop particles that drift below the bathymetry\n", + "\n", + "Because SCHISM's LSC2 grid keeps below-seabed levels as `NaN`, a particle whose (fixed) depth ends up\n", + "deeper than the local water column will sample `NaN`. Rather than feed it a fabricated velocity, we let\n", + "Parcels flag it: sampling `NaN` sets the particle state to `ErrorInterpolation`, and drifting outside\n", + "the mesh sets `ErrorOutOfBounds`. We add a small kernel that runs after advection, records that the\n", + "particle went out of bounds, and sets its state to `Delete` so it stops being advected (its trajectory\n", + "up to that point is kept).\n", + "\n", + "```{note}\n", + "Most of the Lake Ontario mesh is shallow nearshore water with only a couple of valid vertical levels;\n", + "the deep central basin has the full set. We therefore release particles in the deep basin (nodes with\n", + "many valid levels) so they start within the resolved water column.\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "OUT_OF_BOUNDS_STATES = [\n", + " StatusCode.ErrorInterpolation, # sampled NaN (below the local seabed)\n", + " StatusCode.ErrorOutOfBounds, # left the horizontal mesh / below the grid\n", + " StatusCode.ErrorThroughSurface, # above the surface\n", + "]\n", + "\n", + "SchismParticle = parcels.Particle.add_variable(\n", + " parcels.Variable(\"out_of_bounds\", dtype=np.int32, initial=0)\n", + ")\n", + "\n", + "\n", + "def StopBelowBed(particles, fieldset):\n", + " \"\"\"Flag out-of-bounds particles and stop advecting them.\"\"\"\n", + " oob = np.isin(particles.state, OUT_OF_BOUNDS_STATES)\n", + " particles.out_of_bounds = np.where(oob, 1, particles.out_of_bounds)\n", + " particles.state = np.where(oob, StatusCode.Delete, particles.state)" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Release particles and advect" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "# number of valid (non-NaN) vertical levels at each node\n", + "valid_levels = np.isfinite(U.isel(time=0).values).sum(axis=1)\n", + "\n", + "# release on face centroids in the deep basin (all 3 nodes well-resolved vertically)\n", + "cx = node_x[face_nodes].mean(axis=1)\n", + "cy = node_y[face_nodes].mean(axis=1)\n", + "deep = (valid_levels[face_nodes] >= 15).all(axis=1)\n", + "idx = np.where(deep)[0]\n", + "idx = idx[:: max(1, idx.size // 150)][:150]\n", + "\n", + "lon, lat = cx[idx], cy[idx]\n", + "z = np.full(lon.size, 2.0) # release at 2 m depth\n", + "print(f\"releasing {lon.size} particles at z = 2 m in the deep basin\")\n", + "\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset, pclass=SchismParticle, lon=lon, lat=lat, z=z\n", + ")\n", + "output_file = parcels.ParticleFile(\n", + " \"output-schism.parquet\", outputdt=np.timedelta64(30, \"m\")\n", + ")\n", + "\n", + "pset.execute(\n", + " [parcels.kernels.AdvectionRK4, StopBelowBed],\n", + " runtime=np.timedelta64(5, \"h\"), # the dataset spans 5 hours\n", + " dt=np.timedelta64(5, \"m\"),\n", + " output_file=output_file,\n", + " verbose_progress=False,\n", + ")\n", + "print(f\"{len(pset.lon)} of {lon.size} particles still active at the end of the run\")" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "## Plot the velocity field and trajectories\n", + "\n", + "We plot the surface speed across the triangular mesh (in projected kilometres), the particle release\n", + "points, and their trajectories coloured by time since release. The lake currents are slow (~0.1 m/s), so\n", + "over the 5-hour window most particles move only a kilometre or two, while those caught in the faster jet\n", + "near the north-eastern outflow (towards the St. Lawrence) travel noticeably further." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "df = parcels.read_particlefile(\"output-schism.parquet\")\n", + "\n", + "triang = mtri.Triangulation(node_x / 1e3, node_y / 1e3, triangles=face_nodes)\n", + "surf_speed = np.hypot(\n", + " np.asarray(uxds[\"U\"].isel(time=0, zf=0)), np.asarray(uxds[\"V\"].isel(time=0, zf=0))\n", + ")\n", + "speed_face = np.nanmean(surf_speed[face_nodes], axis=1)\n", + "\n", + "fig, ax = plt.subplots(figsize=(11, 7))\n", + "tpc = ax.tripcolor(\n", + " triang,\n", + " facecolors=np.nan_to_num(speed_face),\n", + " shading=\"flat\",\n", + " cmap=\"Blues\",\n", + " vmax=np.nanpercentile(speed_face, 98),\n", + ")\n", + "fig.colorbar(tpc, ax=ax, label=\"surface speed [m/s]\", shrink=0.8)\n", + "\n", + "for traj in df.sort(\"time\").partition_by(\"particle_id\"):\n", + " ax.plot(\n", + " np.array(traj[\"lon\"]) / 1e3,\n", + " np.array(traj[\"lat\"]) / 1e3,\n", + " color=\"0.3\",\n", + " lw=0.6,\n", + " alpha=0.7,\n", + " zorder=2,\n", + " )\n", + "ax.scatter(\n", + " lon / 1e3,\n", + " lat / 1e3,\n", + " facecolors=\"none\",\n", + " edgecolors=\"k\",\n", + " s=20,\n", + " zorder=3,\n", + " label=\"release\",\n", + ")\n", + "\n", + "elapsed_h = (df[\"time\"] - df[\"time\"].min()).dt.total_seconds() / 3600\n", + "sc = ax.scatter(\n", + " np.array(df[\"lon\"]) / 1e3,\n", + " np.array(df[\"lat\"]) / 1e3,\n", + " c=elapsed_h,\n", + " s=4,\n", + " cmap=\"viridis\",\n", + " zorder=3,\n", + ")\n", + "fig.colorbar(sc, ax=ax, label=\"time since release [h]\", shrink=0.8)\n", + "\n", + "ax.set_xlabel(\"projected x [km]\")\n", + "ax.set_ylabel(\"projected y [km]\")\n", + "ax.set_title(\"SCHISM Lake Ontario surface currents with particle trajectories\")\n", + "ax.set_aspect(\"equal\")\n", + "ax.legend(loc=\"upper left\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "The particles drift with the SCHISM surface currents, and any that wander over shallow water where\n", + "their release depth falls below the seabed are flagged (`out_of_bounds == 1`) and stop.\n", + "\n", + "From here, the rest of Parcels works exactly as on structured grids. To go further with SCHISM data:\n", + "\n", + "* Keep in mind that the vertical is approximate: because Parcels uses a single 1D vertical column,\n", + " this tutorial is most meaningful for near-surface and horizontal transport. Faithful full-depth 3D\n", + " transport on an LSC2 grid would require Parcels to support a laterally varying vertical coordinate.\n", + "* Add the vertical velocity as a `W` field and use `AdvectionRK4_3D` for three-dimensional transport\n", + " (still subject to the single-column vertical approximation above).\n", + "* See the [interpolation tutorial](./tutorial_interpolation.ipynb) for the available `Ux*` interpolators." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index 9e22f1236..70cbbea8f 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -29,6 +29,7 @@ examples/tutorial_nemo.ipynb examples/tutorial_croco_3D.ipynb examples/tutorial_mitgcm.ipynb examples/tutorial_fesom.ipynb +examples/tutorial_schism.ipynb examples/tutorial_velocityconversion.ipynb examples/tutorial_nestedgrids.ipynb examples/tutorial_manipulating_field_data.ipynb diff --git a/src/parcels/_datasets/remote.py b/src/parcels/_datasets/remote.py index eb2206591..e176a1fb2 100644 --- a/src/parcels/_datasets/remote.py +++ b/src/parcels/_datasets/remote.py @@ -63,6 +63,11 @@ def _get_data_home() -> Path: "data/FESOM_periodic_channel/v.fesom_channel.nc", "data/FESOM_periodic_channel/w.fesom_channel.nc", ] + + [ + "data/SCHISM_LakeOntario/out2d.schism_lake_ontario.nc", + "data/SCHISM_LakeOntario/horizontalVelX.schism_lake_ontario.nc", + "data/SCHISM_LakeOntario/horizontalVelY.schism_lake_ontario.nc", + ] + [ "data/NemoCurvilinear_data/U_purely_zonal-ORCA025_grid_U.nc4", "data/NemoCurvilinear_data/V_purely_zonal-ORCA025_grid_V.nc4", @@ -222,6 +227,9 @@ class _Purpose(enum.Enum): ("FESOM_periodic_channel/u.fesom_channel", (_V3Dataset(_ODIE,"data/FESOM_periodic_channel/u.fesom_channel.nc"), _Purpose.TUTORIAL)), ("FESOM_periodic_channel/v.fesom_channel", (_V3Dataset(_ODIE,"data/FESOM_periodic_channel/v.fesom_channel.nc"), _Purpose.TUTORIAL)), ("FESOM_periodic_channel/w.fesom_channel", (_V3Dataset(_ODIE,"data/FESOM_periodic_channel/w.fesom_channel.nc"), _Purpose.TUTORIAL)), + ("SCHISM_LakeOntario/out2d", (_V3Dataset(_ODIE,"data/SCHISM_LakeOntario/out2d.schism_lake_ontario.nc"), _Purpose.TUTORIAL)), + ("SCHISM_LakeOntario/horizontalVelX", (_V3Dataset(_ODIE,"data/SCHISM_LakeOntario/horizontalVelX.schism_lake_ontario.nc"), _Purpose.TUTORIAL)), + ("SCHISM_LakeOntario/horizontalVelY", (_V3Dataset(_ODIE,"data/SCHISM_LakeOntario/horizontalVelY.schism_lake_ontario.nc"), _Purpose.TUTORIAL)), ("NemoCurvilinear_data_zonal/U", (_V3Dataset(_ODIE,"data/NemoCurvilinear_data/U_purely_zonal-ORCA025_grid_U.nc4"), _Purpose.TUTORIAL)), ("NemoCurvilinear_data_zonal/V", (_V3Dataset(_ODIE,"data/NemoCurvilinear_data/V_purely_zonal-ORCA025_grid_V.nc4"), _Purpose.TUTORIAL)), ("NemoCurvilinear_data_zonal/mesh_mask", (_V3Dataset(_ODIE,"data/NemoCurvilinear_data/mesh_mask.nc4", _preprocess_drop_time_from_mesh2), _Purpose.TUTORIAL)),