From cc6c18c9221063d9467f6a72f6ee2537edafb01d Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Fri, 12 Jun 2026 11:15:10 -0400 Subject: [PATCH 1/2] Add tutorial for SCHISM model for lake ontario. This is a minimum working example that illustrates primarily how to do 2D (lateral) particle advection with SCHISM model output. It walks through handling issues with uxarray, which assumes that the input grid is in units of degrees for latitude and longitude. The SCHISM output provided here uses x,y in cartesian coordinates in units of meters. Since uxarray automatically wraps longitude to be between [0,360] we have to overwrite the `node_lon` values stored in the uxgrid object. We show how to rename the SCHISM velocity field components, lateral node coordinate, and vertical grid names to match what Parcels expects to see. Note that we intentionally create a ficticious vertical grid in this example; Parcels currently does not support the LSC2 vertical grid that is used in SCHISM and this is documented clearly in the tutorial. --- .../user_guide/examples/tutorial_schism.ipynb | 453 ++++++++++++++++++ docs/user_guide/index.md | 1 + src/parcels/_datasets/remote.py | 8 + 3 files changed, 462 insertions(+) create mode 100644 docs/user_guide/examples/tutorial_schism.ipynb diff --git a/docs/user_guide/examples/tutorial_schism.ipynb b/docs/user_guide/examples/tutorial_schism.ipynb new file mode 100644 index 000000000..0416a4783 --- /dev/null +++ b/docs/user_guide/examples/tutorial_schism.ipynb @@ -0,0 +1,453 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "99ff9ae6", + "metadata": {}, + "source": [ + "# 🖥️ SCHISM tutorial" + ] + }, + { + "cell_type": "markdown", + "id": "0b99ba4b", + "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": "69cd9463", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:08:42.638984Z", + "iopub.status.busy": "2026-06-12T15:08:42.638746Z", + "iopub.status.idle": "2026-06-12T15:09:11.820289Z", + "shell.execute_reply": "2026-06-12T15:09:11.819467Z" + } + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import matplotlib.tri as mtri\n", + "import numpy as np\n", + "import xarray as xr\n", + "import uxarray as ux\n", + "\n", + "import parcels\n", + "import parcels.tutorial\n", + "from parcels._core.statuscodes import StatusCode" + ] + }, + { + "cell_type": "markdown", + "id": "91dc5d03", + "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": "f071e80a", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:09:11.823179Z", + "iopub.status.busy": "2026-06-12T15:09:11.822565Z", + "iopub.status.idle": "2026-06-12T15:09:12.645338Z", + "shell.execute_reply": "2026-06-12T15:09:12.644559Z" + } + }, + "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": "407360fa", + "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": "353da75e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:09:12.647865Z", + "iopub.status.busy": "2026-06-12T15:09:12.647612Z", + "iopub.status.idle": "2026-06-12T15:09:12.855683Z", + "shell.execute_reply": "2026-06-12T15:09:12.854879Z" + } + }, + "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 = grid_ds[\"SCHISM_hgrid_face_nodes\"].values[:, :3].astype(\"int64\") - 1 # 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(f\"n_node={uxgrid.n_node}, n_face={uxgrid.n_face}, n_max_face_nodes={uxgrid.n_max_face_nodes}\")\n", + "print(f\"x range: {node_x.min():.0f} .. {node_x.max():.0f} m\")" + ] + }, + { + "cell_type": "markdown", + "id": "0837cd44", + "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": "75a36eb6", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:09:12.858273Z", + "iopub.status.busy": "2026-06-12T15:09:12.858018Z", + "iopub.status.idle": "2026-06-12T15:09:15.109306Z", + "shell.execute_reply": "2026-06-12T15:09:15.108471Z" + } + }, + "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 = u.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename).load() # 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", + " {\"U\": U.transpose(\"time\", \"zf\", \"n_node\"), \"V\": V.transpose(\"time\", \"zf\", \"n_node\")},\n", + " coords={\"zf\": (\"zf\", zf), \"zc\": (\"zc\", zc)},\n", + " ),\n", + " uxgrid=uxgrid,\n", + ")\n", + "uxds" + ] + }, + { + "cell_type": "markdown", + "id": "7b9ee5f3", + "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": "d3a8bee0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:09:15.111573Z", + "iopub.status.busy": "2026-06-12T15:09:15.111336Z", + "iopub.status.idle": "2026-06-12T15:09:15.119169Z", + "shell.execute_reply": "2026-06-12T15:09:15.118339Z" + } + }, + "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(f\"{name:>4s} -> {type(field).__name__:<11s} interp={interp.__name__ if interp else '-'}\")\n", + "print(\"time interval:\", fieldset.time_interval)" + ] + }, + { + "cell_type": "markdown", + "id": "897796d5", + "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": "9eaf228a", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:09:15.121349Z", + "iopub.status.busy": "2026-06-12T15:09:15.121159Z", + "iopub.status.idle": "2026-06-12T15:09:15.124946Z", + "shell.execute_reply": "2026-06-12T15:09:15.124155Z" + } + }, + "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": "b2275fb8", + "metadata": {}, + "source": [ + "## Release particles and advect" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75820fb3", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:09:15.126973Z", + "iopub.status.busy": "2026-06-12T15:09:15.126794Z", + "iopub.status.idle": "2026-06-12T15:09:16.362213Z", + "shell.execute_reply": "2026-06-12T15:09:16.361455Z" + } + }, + "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(fieldset=fieldset, pclass=SchismParticle, lon=lon, lat=lat, z=z)\n", + "output_file = parcels.ParticleFile(\"output-schism.parquet\", outputdt=np.timedelta64(30, \"m\"))\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": "f4201357", + "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": "55825c33", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-12T15:09:16.364368Z", + "iopub.status.busy": "2026-06-12T15:09:16.364163Z", + "iopub.status.idle": "2026-06-12T15:09:22.680573Z", + "shell.execute_reply": "2026-06-12T15:09:22.679862Z" + } + }, + "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, facecolors=np.nan_to_num(speed_face), shading=\"flat\",\n", + " cmap=\"Blues\", 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(np.array(traj[\"lon\"]) / 1e3, np.array(traj[\"lat\"]) / 1e3, color=\"0.3\", lw=0.6, alpha=0.7, zorder=2)\n", + "ax.scatter(lon / 1e3, lat / 1e3, facecolors=\"none\", edgecolors=\"k\", s=20, zorder=3, label=\"release\")\n", + "\n", + "elapsed_h = (df[\"time\"] - df[\"time\"].min()).dt.total_seconds() / 3600\n", + "sc = ax.scatter(\n", + " np.array(df[\"lon\"]) / 1e3, np.array(df[\"lat\"]) / 1e3, c=elapsed_h, s=4, cmap=\"viridis\", 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": "4d3c1d62", + "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)), From b027ed1e551c03ed72d399d5a99045bb574dffa5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 12 Jun 2026 15:25:12 +0000 Subject: [PATCH 2/2] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../user_guide/examples/tutorial_schism.ipynb | 176 ++++++++---------- 1 file changed, 79 insertions(+), 97 deletions(-) diff --git a/docs/user_guide/examples/tutorial_schism.ipynb b/docs/user_guide/examples/tutorial_schism.ipynb index 0416a4783..432c98e7d 100644 --- a/docs/user_guide/examples/tutorial_schism.ipynb +++ b/docs/user_guide/examples/tutorial_schism.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "99ff9ae6", + "id": "0", "metadata": {}, "source": [ "# 🖥️ SCHISM tutorial" @@ -10,7 +10,7 @@ }, { "cell_type": "markdown", - "id": "0b99ba4b", + "id": "1", "metadata": {}, "source": [ "Parcels v4 supports unstructured-grid model output via [uxarray](https://uxarray.readthedocs.io/).\n", @@ -38,22 +38,15 @@ { "cell_type": "code", "execution_count": null, - "id": "69cd9463", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:08:42.638984Z", - "iopub.status.busy": "2026-06-12T15:08:42.638746Z", - "iopub.status.idle": "2026-06-12T15:09:11.820289Z", - "shell.execute_reply": "2026-06-12T15:09:11.819467Z" - } - }, + "id": "2", + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib.tri as mtri\n", "import numpy as np\n", - "import xarray as xr\n", "import uxarray as ux\n", + "import xarray as xr\n", "\n", "import parcels\n", "import parcels.tutorial\n", @@ -62,7 +55,7 @@ }, { "cell_type": "markdown", - "id": "91dc5d03", + "id": "3", "metadata": {}, "source": [ "## Get the SCHISM tutorial dataset\n", @@ -82,15 +75,8 @@ { "cell_type": "code", "execution_count": null, - "id": "f071e80a", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:09:11.823179Z", - "iopub.status.busy": "2026-06-12T15:09:11.822565Z", - "iopub.status.idle": "2026-06-12T15:09:12.645338Z", - "shell.execute_reply": "2026-06-12T15:09:12.644559Z" - } - }, + "id": "4", + "metadata": {}, "outputs": [], "source": [ "grid_ds = parcels.tutorial.open_dataset(\"SCHISM_LakeOntario/out2d\")\n", @@ -101,7 +87,7 @@ }, { "cell_type": "markdown", - "id": "407360fa", + "id": "5", "metadata": {}, "source": [ "## Build the horizontal mesh\n", @@ -124,20 +110,15 @@ { "cell_type": "code", "execution_count": null, - "id": "353da75e", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:09:12.647865Z", - "iopub.status.busy": "2026-06-12T15:09:12.647612Z", - "iopub.status.idle": "2026-06-12T15:09:12.855683Z", - "shell.execute_reply": "2026-06-12T15:09:12.854879Z" - } - }, + "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 = grid_ds[\"SCHISM_hgrid_face_nodes\"].values[:, :3].astype(\"int64\") - 1 # all-triangular, 0-based\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", @@ -146,13 +127,15 @@ "uxgrid.node_lon.values[:] = node_x\n", "uxgrid.node_lat.values[:] = node_y\n", "\n", - "print(f\"n_node={uxgrid.n_node}, n_face={uxgrid.n_face}, n_max_face_nodes={uxgrid.n_max_face_nodes}\")\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": "0837cd44", + "id": "7", "metadata": {}, "source": [ "## Assemble the velocity fields\n", @@ -181,15 +164,8 @@ { "cell_type": "code", "execution_count": null, - "id": "75a36eb6", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:09:12.858273Z", - "iopub.status.busy": "2026-06-12T15:09:12.858018Z", - "iopub.status.idle": "2026-06-12T15:09:15.109306Z", - "shell.execute_reply": "2026-06-12T15:09:15.108471Z" - } - }, + "id": "8", + "metadata": {}, "outputs": [], "source": [ "nlev = u.sizes[\"nSCHISM_vgrid_layers\"]\n", @@ -199,12 +175,17 @@ "zc = 0.5 * (zf[1:] + zf[:-1])\n", "\n", "rename = {\"nSCHISM_vgrid_layers\": \"zf\", \"nSCHISM_hgrid_node\": \"n_node\"}\n", - "U = u.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename).load() # reverse to surface-first\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", - " {\"U\": U.transpose(\"time\", \"zf\", \"n_node\"), \"V\": V.transpose(\"time\", \"zf\", \"n_node\")},\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", @@ -214,7 +195,7 @@ }, { "cell_type": "markdown", - "id": "7b9ee5f3", + "id": "9", "metadata": {}, "source": [ "## Build the `FieldSet`\n", @@ -229,28 +210,23 @@ { "cell_type": "code", "execution_count": null, - "id": "d3a8bee0", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:09:15.111573Z", - "iopub.status.busy": "2026-06-12T15:09:15.111336Z", - "iopub.status.idle": "2026-06-12T15:09:15.119169Z", - "shell.execute_reply": "2026-06-12T15:09:15.118339Z" - } - }, + "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(f\"{name:>4s} -> {type(field).__name__:<11s} interp={interp.__name__ if interp else '-'}\")\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": "897796d5", + "id": "11", "metadata": {}, "source": [ "## Stop particles that drift below the bathymetry\n", @@ -272,20 +248,13 @@ { "cell_type": "code", "execution_count": null, - "id": "9eaf228a", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:09:15.121349Z", - "iopub.status.busy": "2026-06-12T15:09:15.121159Z", - "iopub.status.idle": "2026-06-12T15:09:15.124946Z", - "shell.execute_reply": "2026-06-12T15:09:15.124155Z" - } - }, + "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.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", @@ -303,7 +272,7 @@ }, { "cell_type": "markdown", - "id": "b2275fb8", + "id": "13", "metadata": {}, "source": [ "## Release particles and advect" @@ -312,15 +281,8 @@ { "cell_type": "code", "execution_count": null, - "id": "75820fb3", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:09:15.126973Z", - "iopub.status.busy": "2026-06-12T15:09:15.126794Z", - "iopub.status.idle": "2026-06-12T15:09:16.362213Z", - "shell.execute_reply": "2026-06-12T15:09:16.361455Z" - } - }, + "id": "14", + "metadata": {}, "outputs": [], "source": [ "# number of valid (non-NaN) vertical levels at each node\n", @@ -337,8 +299,12 @@ "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(fieldset=fieldset, pclass=SchismParticle, lon=lon, lat=lat, z=z)\n", - "output_file = parcels.ParticleFile(\"output-schism.parquet\", outputdt=np.timedelta64(30, \"m\"))\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", @@ -352,7 +318,7 @@ }, { "cell_type": "markdown", - "id": "f4201357", + "id": "15", "metadata": {}, "source": [ "## Plot the velocity field and trajectories\n", @@ -366,15 +332,8 @@ { "cell_type": "code", "execution_count": null, - "id": "55825c33", - "metadata": { - "execution": { - "iopub.execute_input": "2026-06-12T15:09:16.364368Z", - "iopub.status.busy": "2026-06-12T15:09:16.364163Z", - "iopub.status.idle": "2026-06-12T15:09:22.680573Z", - "shell.execute_reply": "2026-06-12T15:09:22.679862Z" - } - }, + "id": "16", + "metadata": {}, "outputs": [], "source": [ "df = parcels.read_particlefile(\"output-schism.parquet\")\n", @@ -387,18 +346,41 @@ "\n", "fig, ax = plt.subplots(figsize=(11, 7))\n", "tpc = ax.tripcolor(\n", - " triang, facecolors=np.nan_to_num(speed_face), shading=\"flat\",\n", - " cmap=\"Blues\", vmax=np.nanpercentile(speed_face, 98),\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(np.array(traj[\"lon\"]) / 1e3, np.array(traj[\"lat\"]) / 1e3, color=\"0.3\", lw=0.6, alpha=0.7, zorder=2)\n", - "ax.scatter(lon / 1e3, lat / 1e3, facecolors=\"none\", edgecolors=\"k\", s=20, zorder=3, label=\"release\")\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, np.array(df[\"lat\"]) / 1e3, c=elapsed_h, s=4, cmap=\"viridis\", zorder=3\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", @@ -412,7 +394,7 @@ }, { "cell_type": "markdown", - "id": "4d3c1d62", + "id": "17", "metadata": {}, "source": [ "The particles drift with the SCHISM surface currents, and any that wander over shallow water where\n",