diff --git a/docs/paper-v4/.gitkeep b/docs/paper-v4/.gitkeep new file mode 100644 index 000000000..e69de29bb diff --git a/docs/paper-v4/data-ingestion.png b/docs/paper-v4/data-ingestion.png new file mode 100644 index 000000000..a3671edef Binary files /dev/null and b/docs/paper-v4/data-ingestion.png differ diff --git a/docs/paper-v4/paper.bib b/docs/paper-v4/paper.bib new file mode 100644 index 000000000..245299f89 --- /dev/null +++ b/docs/paper-v4/paper.bib @@ -0,0 +1,58 @@ +@article{Lange2017, + author = {Lange, Michael and {van Sebille}, Erik}, + doi = {10.5194/gmd-10-4175-2017}, + issn = {1991-9603}, + journal = {Geoscientific Model Development}, + month = {November}, + number = {11}, + pages = {4175–4186}, + publisher = {Copernicus GmbH}, + title = {{Parcels v0.9: prototyping a Lagrangian ocean analysis framework for the petascale age}}, + url = {http://dx.doi.org/10.5194/gmd-10-4175-2017}, + volume = {10}, + year = {2017} +} + +@article{Delandmeter2019, + author = {Delandmeter, Philippe and {van Sebille}, Erik}, + doi = {10.5194/gmd-12-3571-2019}, + issn = {1991-9603}, + journal = {Geoscientific Model Development}, + month = {August}, + number = {8}, + pages = {3571–3584}, + publisher = {Copernicus GmbH}, + title = {{The Parcels v2.0 Lagrangian framework: new field interpolation schemes}}, + url = {http://dx.doi.org/10.5194/gmd-12-3571-2019}, + volume = {12}, + year = {2019} +} + +@article{Hoyer2017, + title = "xarray: {N-D} labeled Arrays and Datasets in Python", + author = "Hoyer, Stephan and Hamman, Joe", + abstract = "xarray is an open source project and Python package that + provides a toolkit and data structures for N-dimensional labeled + arrays. Our approach combines an application programing + interface (API) inspired by pandas with the Common Data Model + for self-described scientific data. Key features of the xarray + package include label-based indexing and arithmetic, + interoperability with the core scientific Python packages (e.g., + pandas, NumPy, Matplotlib), out-of-core computation on datasets + that don't fit into memory, a wide range of serialization and + input/output (I/O) options, and advanced multi-dimensional data + manipulation tools such as group-by and resampling. xarray, as a + data model and analytics toolkit, has been widely adopted in the + geoscience community but is also used more broadly for + multi-dimensional data analysis in physics, machine learning and + finance.", + journal = "J. Open Res. Softw.", + publisher = "Ubiquity Press, Ltd.", + volume = 5, + number = 1, + pages = "10", + month = apr, + year = 2017, + copyright = "http://creativecommons.org/licenses/by/4.0", + language = "en" +} diff --git a/docs/paper-v4/paper.md b/docs/paper-v4/paper.md new file mode 100644 index 000000000..7b3f144df --- /dev/null +++ b/docs/paper-v4/paper.md @@ -0,0 +1,103 @@ +--- +title: "Parcels v4: An Xarray-aligned, flexible lagrangian simulation framework for geoscience" +tags: + - Python + - Lagrangian modelling + - Unstructured grids + - Xarray + - Pangeo +authors: + - name: Nick Hodgskin + orcid: 0009-0003-0778-3183 + corresponding: true + affiliation: 1 + - name: Joe Schoonover + orcid: 0000-0001-5650-7095 + affiliation: 2 + - name: Erik van Sebille + orcid: 0000-0003-2041-0704 + affiliation: 1 + +affiliations: + - name: Institute for Marine and Atmospheric Research, Utrecht University, the Netherlands + index: 1 + - name: Fluid Numerics, Hickory, NC, USA + index: 2 +date: 17 August 2026 +bibliography: paper.bib +--- + +# Summary + +Parcels [@Lange2017; @Delandmeter2019] is a highly customisable Lagrangian simulation framework. +Version 4 of the software is a major update which overhauls the software internals to natively leverage Xarray [@Delandmeter2019] dataset objects. +This makes Parcels compatable with many new data formats (e.g., Zarr, Icechunk) and execution modes (e.g., streaming data from cloud buckets or other data providers - such as the [Copernicus Marine Data Store](https://marine.copernicus.eu/)). +With this update Parcels also adds several new features, including support for unstructured grid datasets (enabling simulations on (combinations of) different grid geometries), support for custom interpolators (surfacing to scientists even more control over the numerics of their simulation), and trajectory output in Parquet format. + +# Statement of need + +The first release of Parcels was in July 2017. +Since then, several factors have changed the climate modelling space. +A significant portion of the geospatial science community has shifted from creating their own scripts for manipulating NetCDF data, to using open source software - namely Xarray - for working with multidimensional climate data. +Major benefits of this software include (a) providing a single in-memory, metadata-rich representation of a full NetCDF-like dataset, (b) its flexibility to manipulate datasets in other data-formats (e.g., Zarr, HDF5), and (c) its abstractions allowing for easy data ingestion from network data sources. + +Writing software that natively works from Xarray datasets provides a powerful abstraction layer allowing downstream developers to create software that works across data formats and data ingestion paradigms. +This is particularly important as climate datasets are continually increasing in resolution and size, preventing local file based storage, and datasets are increasingly being stored in modern data formats such as Zarr. + +Another interesting aspect is that climate modellers are increasingly providing model output on different grid geometries which have attractive features compared to conventional structured grids. +Running natively on these grid geometries, which can be represented as unstructured grids, without re-interpolation allows researchers to run simulations that capture the details from the original dataset. + +Finally, some naming conventions and abstractions in the Parcels codebase were previously oriented towards oceanographers. +These domain specific items have been removed in this version, and our documentation has been adapted so that the applicability of Parcels to other (geo)scientific domains, such as atmospheric or cryospheric particle tracking, is more apparent. + +# State of the field + +> TODO: A description of how this software compares to other commonly-used packages in the research area. If related tools exist, provide a clear “build vs. contribute” justification explaining your unique scholarly contribution and why existing alternatives are insufficient. + +# Software design + +When working with output from a single circulation model, version 4 of Parcels assumes data (i.e., the field data and mesh data) is all contained in a single Xarray dataset object. +This object is either opened directly from disk/a data store, or is constructed by the user from the component Field and mesh files with help from a "converter" function. +These converters also attach relevant CF-convention and grid geometry metadata (SGRID metadata for structured data, UGRID metadata for unstructured data) to the dataset object, allowing the internals of Parcels to assume a certain dataset structure and metadata richness. +Below is a diagram illustrating the code path when working from model data to a fully constructed FieldSet. +If a user wants to run a simulation with fields from different models, they load each model data into its own FieldSet and then combine the FieldSets together into a single FieldSet. + +![](./data-ingestion.png) + +A strength of earlier Parcels versions has been the ability for users to write "kernels" which define particle movement over the course of a simulation. +This added flexibility has enabled users to model a wide range of physical phenomena. +Parcels version 4 adds custom interpolators, allowing users to also have control on how field data is interpolated at particle positions. +Users can use a range of pre-packaged interpolators, or write their own, and set them on a field-by-field basis overriding the default linear interpolators. + +Parcels version 4 also changes the output format from Zarr to Parquet, aligning better with that tabular nature of particle trajectory output. + +# Example use-case + +A key use-case of the Parcels version 4 is the combining of various model data which not only have different sources, but also very different grid geometries. +In this section we present an example simulation focused on the Dutch coast. + +We combine flow data from Deltares’ 3D DCSM-FM model, SWAN Wave model data from Rijkswaterstad, and wind model data from Copernicusmarine. +The 3D DCSM-FM model is an unstructured model, while the other two use structured mesh geometries. +We seed particles along the Dutch coast, and within estuaries facing the North Sea. +We advect the particles for a total of a month (from 2025-11-01 to 2025-12-01) with a timestep of 10 minutes. +The figure below shows the varied types of data and grid geometries that we combine, along with the resulting particle tracks. + +![](./usecase_plot.png) + +# Research impact statement + +Parcels has been cited in over 330 peer reviewed scientific papers so far mostly within the field of oceanography. +This software update expands the reach of Parcels both to more users within oceanography, and in other domains. + +# AI usage disclosure + +Large Language Models (Claude Opus 4.6, ChatGPT ... ) have been used in a guided manner for code generation, documentation, refactoring, and testing. +All LLM written code and documentation has been verified by the authors. +This manuscript was fully written by humans. + +# Acknowledgements + +This work has been funded via Nederlandse Organisatie voor Wetenschappelijk Onderzoek, Exacte en Natuurwetenschappen (VI.C.222.025) as part of the project “Tracing Marine Macroplastics by Unraveling the Ocean’s Multiscale +Transport Processes” and via the Warmworld ELPHE project, and compute on the SURF Research Cloud via project EINF-15719. + +# References diff --git a/docs/paper-v4/usecase_plot.png b/docs/paper-v4/usecase_plot.png new file mode 100644 index 000000000..0f300c75d Binary files /dev/null and b/docs/paper-v4/usecase_plot.png differ diff --git a/docs/paper-v4/usecase_plot.py b/docs/paper-v4/usecase_plot.py new file mode 100644 index 000000000..443c60652 --- /dev/null +++ b/docs/paper-v4/usecase_plot.py @@ -0,0 +1,237 @@ +import cmocean.cm as cmo +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.tri as mtri +import numpy as np +import polars as pl +import xarray as xr + +import parcels + +DIR = "/Users/erik/Desktop/Parcelsv4-paper/" + +# Load Parcels trajecories + +df = parcels.read_particlefile(f"{DIR}/parcels-output-68694.parquet") + +xlim = [df["x"].min(), df["x"].max()] +ylim = [df["y"].min(), df["y"].max()] + +# Build a color map by release time (all particles released together share a color). +release_time_by_pid = df.group_by("particle_id").agg( + pl.col("t").min().alias("release_time") +) +release_times = release_time_by_pid["release_time"].unique().sort().to_list() + +colormap = mpl.colormaps["tab20b"] +release_to_color = { + rt: colormap(i / max(len(release_times) - 1, 1)) + for i, rt in enumerate(release_times) +} +trajectory_to_color = { + row["particle_id"]: release_to_color[row["release_time"]] + for row in release_time_by_pid.iter_rows(named=True) +} + +# Load flow, waves and wind data + +ds_flow = xr.open_dataset(f"{DIR}/dcsm_fm100m_harmonie_202511010000.nc").isel( + time=0, zc=0 +) +zero_faces = (ds_flow["U"] == 0) & (ds_flow["V"] == 0) +ds_flow = ds_flow.isel(n_face=(~zero_faces).values) + +triang = mtri.Triangulation( + ds_flow["Mesh_node_x"].data, + ds_flow["Mesh_node_y"].data, + triangles=ds_flow["tri_face_nodes"].data, +) + +ds_waves = xr.open_dataset(f"{DIR}/swan_kuststrook_harmonie_202511010000.nc").isel( + time=0 +) + +ds_wind = xr.open_dataset(f"{DIR}/copernicusmarine_wind.nc").isel(time=0) + + +# Set up plotting canvas + +speed_xlim = [1, 8] +speed_ylim = [50, 55] +speed_aspect = np.diff(speed_xlim)[0] / np.diff(speed_ylim)[0] +trajectory_aspect = np.diff(xlim)[0] / np.diff(ylim)[0] + +figure_height = 8 +figure_width = figure_height * (trajectory_aspect + speed_aspect / 3) +fig = plt.figure(figsize=(figure_width, figure_height), layout="constrained") +grid = fig.add_gridspec( + 3, + 2, + width_ratios=[speed_aspect, 3 * trajectory_aspect], +) +speed_axes = [fig.add_subplot(grid[row, 0]) for row in range(3)] +particle_ax = fig.add_subplot(grid[:, 1]) + +# Plot particle trajectories and source meshes + +particle_ax.set_xlim(xlim) +particle_ax.set_ylim(ylim) + +for column in range(ds_waves["lon"].values.shape[1]): + label = "Waves mesh" if column == 0 else None + particle_ax.plot( + ds_waves["lon"].values[:, column], + ds_waves["lat"].values[:, column], + color="b", + lw=0.3, + alpha=0.5, + label=label, + ) +for row in range(ds_waves["lon"].values.shape[0]): + particle_ax.plot( + ds_waves["lon"].values[row, :], + ds_waves["lat"].values[row, :], + color="b", + lw=0.3, + alpha=0.5, + ) + +for latitude_index in range(ds_wind["latitude"].values.shape[0]): + label = "Wind mesh" if latitude_index == 0 else None + particle_ax.hlines( + ds_wind["latitude"].values[latitude_index], + xmin=ds_wind["longitude"].values.min(), + xmax=ds_wind["longitude"].values.max(), + color="r", + lw=0.3, + alpha=0.5, + label=label, + zorder=0, + ) +for longitude_index in range(ds_wind["longitude"].values.shape[0]): + particle_ax.vlines( + ds_wind["longitude"].values[longitude_index], + ymin=ds_wind["latitude"].values.min(), + ymax=ds_wind["latitude"].values.max(), + color="r", + lw=0.3, + alpha=0.5, + zorder=1, + ) + +particle_ax.triplot( + triang, + color="k", + lw=0.3, + alpha=0.5, + label="Flow mesh", + zorder=2, +) + +for particle_id in df["particle_id"].unique(): + trajectory = df.filter(pl.col("particle_id") == particle_id) + particle_ax.plot( + trajectory["x"], + trajectory["y"], + color=trajectory_to_color[particle_id], + linewidth=0.6, + alpha=0.3, + ) + particle_ax.plot( + trajectory["x"][-1], + trajectory["y"][-1], + marker="o", + color=trajectory_to_color[particle_id], + markersize=3, + zorder=4, + ) + +legend_handles, legend_labels = particle_ax.get_legend_handles_labels() +particle_ax.legend( + legend_handles[::-1], + legend_labels[::-1], + loc="lower right", +) +legend = particle_ax.get_legend() +for line in legend.get_lines(): + line.set_linewidth(2) + +particle_ax.set_aspect("equal", adjustable="box") +particle_ax.set_title("Parcels trajectories") +particle_ax.set_xlabel("Longitude [°E]") +particle_ax.set_ylabel("Latitude [°N]") + + +# Speed maps +clim = [0, 0.5] +label_box = { + "facecolor": "white", + "edgecolor": "none", + "alpha": 0.8, + "boxstyle": "round,pad=0.3", +} + +# Flow speed +flowspeed = np.hypot(ds_flow["U"], ds_flow["V"]) +flow_map = speed_axes[0].tripcolor( + triang, + facecolors=flowspeed.to_numpy(), + cmap=cmo.speed, + shading="flat", + clim=clim, +) + +# Wave speed +Us = ( + 2 + * np.pi**3 + * ds_waves["wave_height_hm0"] ** 2 + / (ds_waves["wave_period_tm10"] ** 3 * 9.81) +) +ds_waves["Us"] = Us * np.cos(ds_waves["wave_dir_th0"] * np.pi / 180) +ds_waves["Vs"] = Us * np.sin(ds_waves["wave_dir_th0"] * np.pi / 180) +wavespeed = np.hypot(ds_waves["Us"], ds_waves["Vs"]) +wave_map = speed_axes[1].pcolor( + ds_waves["lon"], + ds_waves["lat"], + wavespeed, + cmap=cmo.speed, + clim=clim, +) + +# Wind speed +windspeed = np.hypot(ds_wind["eastward_wind"], ds_wind["northward_wind"]) +wind_map = speed_axes[2].pcolor( + ds_wind["longitude"], + ds_wind["latitude"], + windspeed * 0.01, + cmap=cmo.speed, + clim=clim, +) + +for speed_ax, label in zip(speed_axes, ["Flow", "Waves", "1% of wind"], strict=True): + speed_ax.text( + 0.98, + 0.04, + label, + transform=speed_ax.transAxes, + ha="right", + va="bottom", + bbox=label_box, + ) + speed_ax.set_xlim(speed_xlim) + speed_ax.set_ylim(speed_ylim) + speed_ax.set_aspect("equal", adjustable="box") + +colorbar_ax = speed_axes[2].inset_axes( + [0, -0.34, 1, 0.1], + transform=speed_axes[2].transAxes, + zorder=5, +) +fig.colorbar( + wind_map, + cax=colorbar_ax, + orientation="horizontal", + label="Speed [m s$^{-1}$]", +) +plt.savefig("usecase_plot.png", dpi=300, bbox_inches="tight") diff --git a/docs/paper-v4/usecase_run_script.py b/docs/paper-v4/usecase_run_script.py new file mode 100644 index 000000000..b14066e9e --- /dev/null +++ b/docs/paper-v4/usecase_run_script.py @@ -0,0 +1,223 @@ +import glob +import os + +import copernicusmarine +import numpy as np +import uxarray as ux +import xarray as xr + +import parcels + +DIR = "/storage/shared/oceanparcels/input_data/MatroosWaddenSea/DCSMv7_harmonie" +# %% Open flow files +files = sorted(glob.glob(f"{DIR}/flow/dcsm_fm100m_harmonie_*")) + +ds = xr.open_mfdataset( + files, + combine="nested", + concat_dim="time", + data_vars="minimal", + coords="minimal", + compat="override", + join="override", + parallel=True, + chunks={"time": 1}, +) + +uxgrid = ux.Grid.from_topology( + node_lon=ds["Mesh_node_x"], + node_lat=ds["Mesh_node_y"], + face_node_connectivity=ds["tri_face_nodes"], + fill_value=-1, +) + +uxds = ux.UxDataset( + xr.Dataset( + {"U": ds["U"], "V": ds["V"]}, + coords={ + "zf": ("zf", ds["zf"].data), + "zc": ("zc", ds["zc"].data), + "time": ds["time"], + }, + ), + uxgrid=uxgrid, +) + +fieldset = parcels.FieldSet.from_ugrid_conventions(uxds, mesh="spherical") + +# %% Add Stokes drift to the fieldset +files = sorted(glob.glob(f"{DIR}/waves/swan_kuststrook_harmonie_*.nc")) + +ds = xr.open_mfdataset( + files, + combine="nested", + concat_dim="time", + data_vars="minimal", + coords="minimal", + compat="override", + join="override", + parallel=True, + chunks={"time": 1}, +) +Us = 2 * np.pi**3 * ds["wave_height_hm0"] ** 2 / (ds["wave_period_tm10"] ** 3 * 9.81) +ds["Us"] = Us * np.cos(ds["wave_dir_th0"] * np.pi / 180) +ds["Vs"] = Us * np.sin(ds["wave_dir_th0"] * np.pi / 180) + +ds = ds[["lon", "lat", "Us", "Vs"]].expand_dims(dim={"depth": [0]}) + +for var in ["lon", "lat"]: + ds[var] = ds[var].transpose("col", "row") + +ds["grid"] = xr.DataArray( + 0, + attrs=parcels._sgrid.SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("row", "col"), + node_coordinates=("lon", "lat"), + face_dimensions=( + parcels._sgrid.FaceNodePadding("X", "row", parcels._sgrid.Padding.LOW), + parcels._sgrid.FaceNodePadding("Y", "col", parcels._sgrid.Padding.LOW), + ), + vertical_dimensions=( + parcels._sgrid.FaceNodePadding("Z", "depth", parcels._sgrid.Padding.HIGH), + ), + ).to_attrs(), +) + +fieldset_waves = parcels.FieldSet.from_sgrid_conventions( + ds, vector_fields={"UVStokes": ("Us", "Vs")} +) +fieldset += fieldset_waves + + +# %% Add wind to the fieldset + +startdate = np.datetime64("2025-11-01T00:00:00") +enddate = np.datetime64("2025-12-01T00:00:00") +ds = copernicusmarine.open_dataset( + dataset_id="cmems_obs-wind_glo_phy_my_l4_0.125deg_PT1H", + variables=["eastward_wind", "northward_wind"], + minimum_longitude=1, + maximum_longitude=8, + minimum_latitude=51, + maximum_latitude=55, + start_datetime=np.datetime_as_string(startdate, unit="s"), + end_datetime=np.datetime_as_string(enddate, unit="s"), +) +ds = parcels.convert.copernicusmarine_to_sgrid( + fields={ + "eastward_wind": ds["eastward_wind"], + "northward_wind": ds["northward_wind"], + } +) +ds.load() # Data is mall enough to load into memory +fieldset_wind = parcels.FieldSet.from_sgrid_conventions( + ds, vector_fields={"UVWind": ("eastward_wind", "northward_wind")} +) +fieldset += fieldset_wind + +fieldset = fieldset.to_windowed_arrays() +fieldset.describe() + +# %% Create the simulation +release = "coast" # 'coast' or 'off_shore' + +fieldset.windage = 0.01 # windage factor for the particles + +if release == "off_shore": + lat0 = 52.10 + lon0 = 3.7 +elif release == "coast": + lat0 = 52.020000 + lon0 = 4.097500 + +radii = [100, 200, 400, 800, 1600] # m +n_points = 16 +R = 6371000 # Earth radius (m) + +lat = [lat0] +lon = [lon0] +for r in radii: + theta = np.linspace(0, 2 * np.pi, n_points, endpoint=False) + dlat = np.rad2deg((r / R) * np.cos(theta)) + dlon = np.rad2deg((r / (R * np.cos(np.deg2rad(lat0)))) * np.sin(theta)) + lat.extend(lat0 + dlat) + lon.extend(lon0 + dlon) + +release_dt = np.timedelta64(12, "h") +nrepeat = np.timedelta64(28, "D") // release_dt # number of releases +npart = len(lon) +lon = np.broadcast_to(lon, (nrepeat, npart)) +lat = np.broadcast_to(lat, (nrepeat, npart)) +time_i = fieldset.time_interval.left +time = ( + np.broadcast_to(time_i, (nrepeat, npart)) + + np.arange(0, nrepeat)[:, np.newaxis] * release_dt +) +print( + f"Running {nrepeat} releases of {npart} particles each, for a total of {nrepeat * npart} particles." +) + +MatroosParticle = parcels.Particle.add_variable( + parcels.Variable("outside_stokes", dtype=np.int32, initial=0.0) +) +pset = parcels.ParticleSet(fieldset, pclass=MatroosParticle, x=lon, y=lat, t=time) + +slurm_job_id = os.getenv("SLURM_JOB_ID", "local") +output_name = f"parcels-output-{slurm_job_id}.parquet" + +outputdt = np.timedelta64(30, "m") # output every 30 minutes +output_file = parcels.ParticleFile( + output_name, + outputdt=outputdt, + mode="w", +) + + +def DeleteAnyError(particles, fieldset): + any_error = particles.state >= 50 # This captures all Errors + particles[any_error].state = parcels.StatusCode.Delete + + +def AdvectionRK2(particles, fieldset): # pragma: no cover + """Advection of particles using second-order Runge-Kutta integration.""" + (u1, v1) = fieldset.UV[particles] + (us1, vs1) = fieldset.UVStokes[particles] + (uw1, vw1) = fieldset.UVWind[particles] + x1 = particles.x + (u1 + us1 + uw1 * fieldset.windage) * 0.5 * particles.dt + y1 = particles.y + (v1 + vs1 + vw1 * fieldset.windage) * 0.5 * particles.dt + + (u2, v2) = fieldset.UV[ + particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles + ] + (us2, vs2) = fieldset.UVStokes[ + particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles + ] + (uw2, vw2) = fieldset.UVWind[ + particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles + ] + + # Handle particles that are outside the Stokes drift field + outside_stokes = (us2 == 0) | (vs2 == 0) + us2[outside_stokes] = 0.0 + vs2[outside_stokes] = 0.0 + particles.state[outside_stokes] = parcels.StatusCode.Evaluate + particles.outside_stokes[outside_stokes] = 1 + particles.outside_stokes[~outside_stokes] = 0 + + # set wind to zero for particles that are on land + on_land = (u2 == 0) & (v2 == 0) + uw2[on_land] = 0.0 + vw2[on_land] = 0.0 + + particles.dx += (u2 + us2 + uw2 * fieldset.windage) * particles.dt + particles.dy += (v2 + vs2 + vw2 * fieldset.windage) * particles.dt + + +pset.execute( + [AdvectionRK2, DeleteAnyError], + endtime=fieldset.time_interval.right, + dt=np.timedelta64(10, "m"), + output_file=output_file, +)