From 18bd8b4c7feb5c756c97ca69808d64b2701c4d0b Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 09:28:23 +0100 Subject: [PATCH 01/38] First implementation of fieldset.from_nemo And update of tutorial_nemo_curvilinear --- .../examples/tutorial_nemo_curvilinear.ipynb | 56 +++---- src/parcels/_core/fieldset.py | 144 +++++++++++++++--- 2 files changed, 142 insertions(+), 58 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index 7eb6b0081..57c42b0e1 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb @@ -64,47 +64,23 @@ "outputs": [], "source": [ "data_folder = parcels.download_example_dataset(\"NemoCurvilinear_data\")\n", - "files = data_folder.glob(\"*.nc4\")\n", "ds_fields = xr.open_mfdataset(\n", - " files, combine=\"nested\", data_vars=\"minimal\", coords=\"minimal\", compat=\"override\"\n", + " data_folder.glob(\"*.nc4\"),\n", + " combine=\"nested\",\n", + " data_vars=\"minimal\",\n", + " coords=\"minimal\",\n", + " compat=\"override\",\n", ")\n", - "\n", - "\n", - "# TODO: replace manual fieldset creation with FieldSet.from_nemo() once available\n", "ds_fields = (\n", - " ds_fields.isel(time_counter=0, drop=True)\n", - " .isel(time=0, drop=True)\n", - " .isel(z_a=0, drop=True)\n", - " .rename({\"glamf\": \"lon\", \"gphif\": \"lat\", \"z\": \"depth\"})\n", + " ds_fields.isel(time=0, drop=True).isel(z_a=0, drop=True).isel(z=0, drop=True)\n", ")\n", "\n", - "import xgcm\n", - "\n", - "xgcm_grid = xgcm.Grid(\n", - " ds_fields,\n", - " coords={\n", - " \"X\": {\"left\": \"x\"},\n", - " \"Y\": {\"left\": \"y\"},\n", - " },\n", - " periodic=False,\n", - " autoparse_metadata=False,\n", - ")\n", - "grid = parcels.XGrid(xgcm_grid, mesh=\"spherical\")\n", + "# TODO support 2D spatial fields in NEMO\n", + "ds_fields = ds_fields.expand_dims({\"depth\": [0.0]})\n", + "ds_fields[\"glamf\"] = ds_fields[\"glamf\"].isel(depth=0, drop=True)\n", + "ds_fields[\"gphif\"] = ds_fields[\"gphif\"].isel(depth=0, drop=True)\n", "\n", - "U = parcels.Field(\n", - " \"U\", ds_fields[\"U\"], grid, interp_method=parcels.interpolators.XLinear\n", - ")\n", - "V = parcels.Field(\n", - " \"V\", ds_fields[\"V\"], grid, interp_method=parcels.interpolators.XLinear\n", - ")\n", - "U.units = parcels.GeographicPolar()\n", - "V.units = (\n", - " parcels.GeographicPolar()\n", - ") # U and V need GeographicPolar for C-Grid interpolation to work correctly\n", - "UV = parcels.VectorField(\n", - " \"UV\", U, V, vector_interp_method=parcels.interpolators.CGrid_Velocity\n", - ")\n", - "fieldset = parcels.FieldSet([U, V, UV])" + "fieldset = parcels.FieldSet.from_nemo(ds_fields)" ] }, { @@ -122,7 +98,7 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n", - "pc1 = ds_fields.U.plot(cmap=\"viridis\", ax=ax[0], vmin=0)\n", + "pc1 = fieldset.U.data.plot(cmap=\"viridis\", ax=ax[0], vmin=0)\n", "pc2 = ax[1].pcolormesh(\n", " fieldset.U.grid.lon,\n", " fieldset.U.grid.lat,\n", @@ -152,8 +128,12 @@ "outputs": [], "source": [ "u, v = fieldset.UV.eval(\n", - " np.array([0]), np.array([0]), np.array([20]), np.array([50]), applyConversion=False\n", - ") # do not convert m/s to deg/s\n", + " time=np.array([0]),\n", + " z=np.array([0]),\n", + " y=np.array([40]),\n", + " x=np.array([50]),\n", + " applyConversion=False, # do not convert m/s to deg/s\n", + ")\n", "print(f\"(u, v) = ({u[0]:.3f}, {v[0]:.3f})\")\n", "assert np.isclose(u, 1.0, atol=1e-3)" ] diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index f9bcb795f..4f863fcf0 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -19,7 +19,13 @@ from parcels._core.xgrid import _DEFAULT_XGCM_KWARGS, XGrid from parcels._logger import logger from parcels._typing import Mesh -from parcels.interpolators import UxPiecewiseConstantFace, UxPiecewiseLinearNode, XConstantField, XLinear +from parcels.interpolators import ( + CGrid_Velocity, + UxPiecewiseConstantFace, + UxPiecewiseLinearNode, + XConstantField, + XLinear, +) if TYPE_CHECKING: from parcels._core.basegrid import BaseGrid @@ -205,7 +211,7 @@ def from_copernicusmarine(ds: xr.Dataset): """ ds = ds.copy() - ds = _discover_copernicusmarine_U_and_V(ds) + ds = _discover_U_and_V(ds, _COPERNICUS_MARINE_CF_STANDARD_NAME_FALLBACKS) expected_axes = set("XYZT") # TODO: Update after we have support for 2D spatial fields if missing_axes := (expected_axes - set(ds.cf.axes)): raise ValueError( @@ -215,7 +221,7 @@ def from_copernicusmarine(ds: xr.Dataset): "HINT: Add xarray metadata attribute 'axis' to dimension - e.g., ds['lat'].attrs['axis'] = 'Y'" ) - ds = _rename_coords_copernicusmarine(ds) + ds = _maybe_rename_coords(ds, _COPERNICUS_MARINE_AXIS_VARNAMES) if "W" in ds.data_vars: # Negate W to convert from up positive to down positive (as that's the direction of positive z) ds["W"] -= ds["W"] @@ -239,6 +245,48 @@ def from_copernicusmarine(ds: xr.Dataset): ) return FieldSet.from_sgrid_conventions(ds, mesh="spherical") + def from_nemo(ds: xr.Dataset): + ds = ds.copy() + ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) + ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) + ds = _maybe_rename_coords(ds, _NEMO_AXIS_VARNAMES) + ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_NAMES) + ds = _set_coords(ds, _NEMO_DIMENSION_NAMES) + ds = _set_axis_attrs(ds, _NEMO_AXIS_VARNAMES) + + expected_axes = set("XYZT") # TODO: Update after we have support for 2D spatial fields + if missing_axes := (expected_axes - set(ds.cf.axes)): + raise ValueError( + f"Dataset missing CF compliant metadata for axes " + f"{missing_axes}. Expected 'axis' attribute to be set " + f"on all dimension axes {expected_axes}. " + "HINT: Add xarray metadata attribute 'axis' to dimension - e.g., ds['lat'].attrs['axis'] = 'Y'" + ) + + if "W" in ds.data_vars: + # Negate W to convert from up positive to down positive (as that's the direction of positive z) + ds["W"] -= ds["W"] + if "grid" in ds.cf.cf_roles: + raise ValueError( + "Dataset already has a 'grid' variable (according to cf_roles). Didn't expect there to be grid metadata on copernicusmarine datasets - please open an issue with more information about your dataset." + ) + ds["grid"] = xr.DataArray( + 0, + attrs=sgrid.Grid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("lon", "lat"), + face_dimensions=( + sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), + sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), + ), + vertical_dimensions=(sgrid.DimDimPadding("depth_center", "depth", sgrid.Padding.LOW),), + ).to_attrs(), + ) + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") + fieldset.UV.vector_interp_method = CGrid_Velocity + return fieldset + def from_fesom2(ds: ux.UxDataset): """Create a FieldSet from a FESOM2 uxarray.UxDataset. @@ -404,19 +452,8 @@ def _format_calendar_error_message(field: Field, reference_datetime: TimeLike) - } -def _rename_coords_copernicusmarine(ds): - try: - for axis, [coord] in ds.cf.axes.items(): - ds = ds.rename({coord: _COPERNICUS_MARINE_AXIS_VARNAMES[axis]}) - except ValueError as e: - raise ValueError(f"Multiple coordinates found for Copernicus dataset on axis '{axis}'. Check your data.") from e - return ds - - -def _discover_copernicusmarine_U_and_V(ds: xr.Dataset) -> xr.Dataset: - # Assumes that the dataset has U and V data - - cf_UV_standard_name_fallbacks = [ +_COPERNICUS_MARINE_CF_STANDARD_NAME_FALLBACKS = { + "UV": [ ( "eastward_sea_water_velocity", "northward_sea_water_velocity", @@ -438,11 +475,78 @@ def _discover_copernicusmarine_U_and_V(ds: xr.Dataset) -> xr.Dataset: "eastward_sea_water_velocity_vertical_mean_over_pelagic_layer", "northward_sea_water_velocity_vertical_mean_over_pelagic_layer", ), # GLOBAL_MULTIYEAR_BGC_001_033 - ] - cf_W_standard_name_fallbacks = ["upward_sea_water_velocity", "vertical_sea_water_velocity"] + ], + "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], +} + +_NEMO_CF_STANDARD_NAME_FALLBACKS = { + "UV": [ + ( + "sea_water_x_velocity", + "sea_water_y_velocity", + ), + ], + "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], +} + +_NEMO_DIMENSION_NAMES = ["x", "y", "depth", "time", "lon", "lat"] + +_NEMO_AXIS_VARNAMES = { + "X": "lon", + "Y": "lat", + "Z": "depth", + "T": "time", +} + +_NEMO_VARNAMES_MAPPING = { + "glamf": "lon", + "gphif": "lat", + "time_counter": "time", +} + + +def _maybe_rename_coords(ds, AXIS_VARNAMES): + try: + for axis, [coord] in ds.cf.axes.items(): + ds = ds.rename({coord: AXIS_VARNAMES[axis]}) + except ValueError as e: + raise ValueError(f"Multiple coordinates found for Copernicus dataset on axis '{axis}'. Check your data.") from e + return ds + + +def _maybe_rename_variables(ds, VARNAMES_MAPPING): + rename_dict = {old: new for old, new in VARNAMES_MAPPING.items() if (old in ds.data_vars) or (old in ds.coords)} + if rename_dict: + ds = ds.rename(rename_dict) + return ds + + +def _assign_dims_as_coords(ds, DIMENSION_NAMES): + for axis in DIMENSION_NAMES: + if axis in ds.dims and axis not in ds.coords: + ds = ds.assign_coords({axis: np.arange(ds.sizes[axis])}) + return ds + + +def _set_coords(ds, DIMENSION_NAMES): + for varname in DIMENSION_NAMES: + if varname in ds and varname not in ds.coords: + ds = ds.set_coords([varname]) + return ds + + +def _set_axis_attrs(ds, AXIS_VARNAMES): + for axis, varname in AXIS_VARNAMES.items(): + if varname in ds.coords: + ds[varname].attrs["axis"] = axis + return ds + + +def _discover_U_and_V(ds: xr.Dataset, cf_standard_names_fallbacks) -> xr.Dataset: + # Assumes that the dataset has U and V data if "W" not in ds: - for cf_standard_name_W in cf_W_standard_name_fallbacks: + for cf_standard_name_W in cf_standard_names_fallbacks["W"]: if cf_standard_name_W in ds.cf.standard_names: ds = _ds_rename_using_standard_names(ds, {cf_standard_name_W: "W"}) break @@ -454,7 +558,7 @@ def _discover_copernicusmarine_U_and_V(ds: xr.Dataset) -> xr.Dataset: "Dataset has only one of the two variables 'U' and 'V'. Please rename the appropriate variable in your dataset to have both 'U' and 'V' for Parcels simulation." ) - for cf_standard_name_U, cf_standard_name_V in cf_UV_standard_name_fallbacks: + for cf_standard_name_U, cf_standard_name_V in cf_standard_names_fallbacks["UV"]: if cf_standard_name_U in ds.cf.standard_names: if cf_standard_name_V not in ds.cf.standard_names: raise ValueError( From 57d4a5451b6b12e3add84d70475a18bdba3cf4e1 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 11:23:16 +0100 Subject: [PATCH 02/38] Setting vector_interp_method only when vector fields available --- src/parcels/_core/fieldset.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 4f863fcf0..61394e7fc 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -284,7 +284,10 @@ def from_nemo(ds: xr.Dataset): ).to_attrs(), ) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") - fieldset.UV.vector_interp_method = CGrid_Velocity + if "UV" in fieldset.fields: + fieldset.UV.vector_interp_method = CGrid_Velocity + if "UVW" in fieldset.fields: + fieldset.UVW.vector_interp_method = CGrid_Velocity return fieldset def from_fesom2(ds: ux.UxDataset): @@ -510,7 +513,7 @@ def _maybe_rename_coords(ds, AXIS_VARNAMES): for axis, [coord] in ds.cf.axes.items(): ds = ds.rename({coord: AXIS_VARNAMES[axis]}) except ValueError as e: - raise ValueError(f"Multiple coordinates found for Copernicus dataset on axis '{axis}'. Check your data.") from e + raise ValueError(f"Multiple coordinates found on axis '{axis}'. Check your DataSet.") from e return ds From d47400b1f8680fcfcd53557be6fece036e3a7f43 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 11:37:14 +0100 Subject: [PATCH 03/38] Add error when pset.execute has neither runtime nor endtime --- src/parcels/_core/particleset.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parcels/_core/particleset.py b/src/parcels/_core/particleset.py index 8dea8efaa..839ab0fbd 100644 --- a/src/parcels/_core/particleset.py +++ b/src/parcels/_core/particleset.py @@ -547,6 +547,9 @@ def _get_simulation_start_and_end_times( if runtime is None and time_interval is None: raise ValueError("The runtime must be provided when the time_interval is not defined for a fieldset.") + if runtime is None and endtime is None: + raise ValueError("Either runtime or endtime must be provided.") + if sign_dt == 1: first_release_time = particle_release_times.min() else: From a60c9119a49f156ddecb695ed526ed133cf8df51 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 16:52:58 +0100 Subject: [PATCH 04/38] Adding support for NEMO 3D And updating tutorial --- .../examples/tutorial_nemo_curvilinear.ipynb | 7 +- .../examples_v3/tutorial_nemo_3D.ipynb | 183 ++++++------------ docs/user_guide/index.md | 2 +- src/parcels/_core/fieldset.py | 33 +++- 4 files changed, 92 insertions(+), 133 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index 57c42b0e1..6c980f03c 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb @@ -66,7 +66,6 @@ "data_folder = parcels.download_example_dataset(\"NemoCurvilinear_data\")\n", "ds_fields = xr.open_mfdataset(\n", " data_folder.glob(\"*.nc4\"),\n", - " combine=\"nested\",\n", " data_vars=\"minimal\",\n", " coords=\"minimal\",\n", " compat=\"override\",\n", @@ -76,9 +75,9 @@ ")\n", "\n", "# TODO support 2D spatial fields in NEMO\n", - "ds_fields = ds_fields.expand_dims({\"depth\": [0.0]})\n", - "ds_fields[\"glamf\"] = ds_fields[\"glamf\"].isel(depth=0, drop=True)\n", - "ds_fields[\"gphif\"] = ds_fields[\"gphif\"].isel(depth=0, drop=True)\n", + "ds_fields = ds_fields.expand_dims({\"depthw\": [0.0]})\n", + "ds_fields[\"glamf\"] = ds_fields[\"glamf\"].isel(depthw=0, drop=True)\n", + "ds_fields[\"gphif\"] = ds_fields[\"gphif\"].isel(depthw=0, drop=True)\n", "\n", "fieldset = parcels.FieldSet.from_nemo(ds_fields)" ] diff --git a/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb b/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb index 21af7ad97..805a293a8 100644 --- a/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb +++ b/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb @@ -18,26 +18,13 @@ "More information about C-grid interpolation can be found in [Delandmeter et al., 2019](https://www.geosci-model-dev-discuss.net/gmd-2018-339/).\n", "An example of such a discretisation is the NEMO model, which is one of the models supported in Parcels. A tutorial teaching how to [interpolate 2D data on a NEMO grid](https://docs.oceanparcels.org/en/latest/examples/tutorial_nemo_curvilinear.html) is available within Parcels.\n", "\n", - "Here, we focus on 3D fields. Basically, it is a straightforward extension of the 2D example, but it is very easy to do a mistake in the setup of the vertical discretisation that would affect the interpolation scheme.\n", - "\n", - "## Preliminary comments\n", - "\n", + "```{note}\n", "_How to know if your data is discretised on a C grid?_ The best way is to read the documentation that comes with the data. Alternatively, an easy check is to assess the coordinates of the U, V and W fields: for an A grid, U, V and W are distributed on the same nodes, such that the coordinates are the same. For a C grid, there is a shift of half a cell between the different variables.\n", + "```\n", "\n", - "_What about grid indexing?_ Since the C-grid variables are not located on the same nodes, there is not one obvious way to define the indexing, i.e. where is `u[k,j,i]` compared to `v[k,j,i]` and `w[k,j,i]`. In Parcels, we use the same notation as in NEMO: see [horizontal indexing](https://www.nemo-ocean.eu/doc/img360.png) and [vertical indexing](https://www.nemo-ocean.eu/doc/img362.png).\n", - "It is important that you check if your data is following the same notation. Otherwise, you should re-index your data properly (this can be done within Parcels, there is no need to regenerate new netcdf files).\n", - "\n", - "_What about the accuracy?_ By default in Parcels, particle coordinates (i.e. longitude, latitude and depth) are stored using single-precision `np.float32` numbers. The advantage of this is that it saves some memory resources for the computation. In some applications, especially where particles travel very close to the coast, the single precision accuracy can lead to uncontrolled particle beaching, due to numerical rounding errors. In such case, you may want to double the coordinate precision to `np.float64`. This can be done by adding the parameter `lonlatdepth_dtype=np.float64` to the ParticleSet constructor. Note that for C-grid fieldsets such as in NEMO, the coordinates precision is set by default to `np.float64`.\n", - "\n", - "## How to create a 3D NEMO `dimensions` dictionary?\n", - "\n", - "In the following, we will show how to create the `dimensions` dictionary for 3D NEMO simulations. What you require is a 'mesh_mask' file, which in our case is called `coordinates.nc` but in some other versions of NEMO has a different name. In any case, it will have to contain the variables `glamf`, `gphif` and `depthw`, which are the longitude, latitude and depth of the mesh nodes, respectively. Note that `depthw` is not part of the mesh_mask file, but is in the same file as the w data (`wfiles[0]`).\n", - "\n", - "For the C-grid interpolation in Parcels to work properly, it is important that `U`, `V` and `W` are on the same grid. \n", + "Here, we focus on 3D fields. Basically, it is a straightforward extension of the 2D example, but it is very easy to make a mistake in the setup of the vertical discretisation that would affect the interpolation scheme.\n", "\n", - "All other tracers (e.g. `temperature`, `salinity`) should also be on this same grid. So even though these tracers are computed by NEMO on the T-points, Parcels expects them on the f-points (`glamf`, `gphif` and `depthw`). Parcels then under the hood makes sure the interpolation of these tracers is done correctly.\n", - "\n", - "The code below is an example of how to create a 3D simulation with particles, starting in the mouth of the river Rhine at 1m depth, and advecting them through the North Sea using the `AdvectionRK4_3D`\n" + "For the C-grid interpolation in Parcels to work properly, it is important that `U`, `V` and `W` are on the same grid. All other tracers (e.g. `temperature`, `salinity`) should also be on this same grid. So even though these tracers are computed by NEMO on the T-points, Parcels expects them on the f-points (`glamf`, `gphif` and `depthw`). Parcels then under the hood makes sure the interpolation of these tracers is done correctly." ] }, { @@ -46,69 +33,18 @@ "metadata": {}, "outputs": [], "source": [ - "import warnings\n", - "from datetime import timedelta\n", - "from glob import glob\n", - "\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "import xarray as xr\n", "\n", - "import parcels\n", - "from parcels import FileWarning\n", - "\n", - "# Add a filter for the xarray decoding warning\n", - "warnings.simplefilter(\"ignore\", FileWarning)\n", - "\n", - "example_dataset_folder = parcels.download_example_dataset(\n", - " \"NemoNorthSeaORCA025-N006_data\"\n", - ")\n", - "ufiles = sorted(glob(f\"{example_dataset_folder}/ORCA*U.nc\"))\n", - "vfiles = sorted(glob(f\"{example_dataset_folder}/ORCA*V.nc\"))\n", - "wfiles = sorted(glob(f\"{example_dataset_folder}/ORCA*W.nc\"))\n", - "# tfiles = sorted(glob(f\"{example_dataset_folder}/ORCA*T.nc\")) # Not used in this example\n", - "mesh_mask = f\"{example_dataset_folder}/coordinates.nc\"\n", - "\n", - "filenames = {\n", - " \"U\": {\"lon\": mesh_mask, \"lat\": mesh_mask, \"depth\": wfiles[0], \"data\": ufiles},\n", - " \"V\": {\"lon\": mesh_mask, \"lat\": mesh_mask, \"depth\": wfiles[0], \"data\": vfiles},\n", - " \"W\": {\"lon\": mesh_mask, \"lat\": mesh_mask, \"depth\": wfiles[0], \"data\": wfiles},\n", - " # \"T\": {\"lon\": mesh_mask, \"lat\": mesh_mask, \"depth\": wfiles[0], \"data\": tfiles}, # Not used in this example\n", - "}\n", - "\n", - "variables = {\n", - " \"U\": \"uo\",\n", - " \"V\": \"vo\",\n", - " \"W\": \"wo\",\n", - " # \"T\": \"thetao\", # Not used in this example\n", - "}\n", - "\n", - "# Note that all variables need the same dimensions in a C-Grid\n", - "c_grid_dimensions = {\n", - " \"lon\": \"glamf\",\n", - " \"lat\": \"gphif\",\n", - " \"depth\": \"depthw\",\n", - " \"time\": \"time_counter\",\n", - "}\n", - "dimensions = {\n", - " \"U\": c_grid_dimensions,\n", - " \"V\": c_grid_dimensions,\n", - " \"W\": c_grid_dimensions,\n", - " # \"T\": c_grid_dimensions, # Not used in this example\n", - "}\n", - "\n", - "fieldset = parcels.FieldSet.from_nemo(filenames, variables, dimensions)\n", - "\n", - "pset = parcels.ParticleSet.from_line(\n", - " fieldset=fieldset,\n", - " pclass=parcels.Particle,\n", - " size=10,\n", - " start=(1.9, 52.5),\n", - " finish=(3.4, 51.6),\n", - " depth=1,\n", - ")\n", - "\n", - "kernels = pset.Kernel(parcels.kernels.AdvectionRK4_3D)\n", - "pset.execute(kernels, runtime=timedelta(days=4), dt=timedelta(hours=6))" + "import parcels" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parcels v4 comes with a `from_nemo` method that automatically sets up the correct `Grid` for both 2D and 3D NEMO data. However, you may need to do some data-wrangling on your NEMO data files, such as below:" ] }, { @@ -117,47 +53,65 @@ "metadata": {}, "outputs": [], "source": [ - "depth_level = 8\n", - "print(\n", - " f\"Level[{int(depth_level)}] depth is: \"\n", - " f\"[{fieldset.W.grid.depth[depth_level]:g} \"\n", - " f\"{fieldset.W.grid.depth[depth_level + 1]:g}]\"\n", + "data_folder = parcels.download_example_dataset(\"NemoNorthSeaORCA025-N006_data\")\n", + "\n", + "ds_fields = xr.open_mfdataset(\n", + " data_folder.glob(\"ORCA*.nc\"),\n", + " data_vars=\"minimal\",\n", + " coords=\"minimal\",\n", + " compat=\"override\",\n", ")\n", "\n", - "plt.pcolormesh(\n", - " fieldset.U.grid.lon, fieldset.U.grid.lat, fieldset.U.data[0, depth_level, :, :]\n", + "ds_coord = xr.open_dataset(data_folder / \"coordinates.nc\", decode_times=False).isel(\n", + " time=0, drop=True\n", ")\n", - "plt.show()" + "ds_fields = xr.merge([ds_fields, ds_coord[[\"glamf\", \"gphif\"]]])\n", + "\n", + "fieldset = parcels.FieldSet.from_nemo(ds_fields)" ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "## Adding other fields like cell edges\n", - "\n", - "It is quite straightforward to add other gridded data, on the same curvilinear or any other type of grid, to the fieldset. Because it is good practice to make no changes to a `FieldSet` once a `ParticleSet` has been defined in it, we redefine the fieldset and add the fields with the cell edges from the coordinates file using `FieldSet.add_field()`.\n", - "\n", - "Note that including tracers like `temperature` and `salinity` needs to be done at the f-points, so on the same grid as the velocity fields. See also the section above.\n" + "The code below is an example of how to then create a 3D simulation with particles, starting on a meridional line through the North Sea from the mouth of the river Rhine at 1m depth, and advecting them through the North Sea using the `AdvectionRK2_3D`" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "tags": [ + "hide-output" + ] + }, "outputs": [], "source": [ - "from parcels import Field" + "npart = 10\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " lon=np.linspace(1.9, 3.4, npart),\n", + " lat=np.linspace(65, 51.6, npart),\n", + " z=100 * np.ones(npart),\n", + ")\n", + "\n", + "pfile = parcels.ParticleFile(\n", + " store=\"output_nemo3D.zarr\", outputdt=np.timedelta64(1, \"D\")\n", + ")\n", + "\n", + "pset.execute(\n", + " parcels.kernels.AdvectionRK2_3D,\n", + " endtime=fieldset.time_interval.right,\n", + " dt=np.timedelta64(6, \"h\"),\n", + " output_file=pfile,\n", + ")" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "fieldset = parcels.FieldSet.from_nemo(filenames, variables, dimensions)" + "We can then plot the trajectories on top of the surface U field" ] }, { @@ -166,34 +120,13 @@ "metadata": {}, "outputs": [], "source": [ - "e1u = Field.from_netcdf(\n", - " filenames=mesh_mask,\n", - " variable=\"e1u\",\n", - " dimensions={\"lon\": \"glamu\", \"lat\": \"gphiu\"},\n", - " interp_method=\"nearest\",\n", - ")\n", - "e2u = Field.from_netcdf(\n", - " filenames=mesh_mask,\n", - " variable=\"e2u\",\n", - " dimensions={\"lon\": \"glamu\", \"lat\": \"gphiu\"},\n", - " interp_method=\"nearest\",\n", - ")\n", - "e1v = Field.from_netcdf(\n", - " filenames=mesh_mask,\n", - " variable=\"e1v\",\n", - " dimensions={\"lon\": \"glamv\", \"lat\": \"gphiv\"},\n", - " interp_method=\"nearest\",\n", - ")\n", - "e2v = Field.from_netcdf(\n", - " filenames=mesh_mask,\n", - " variable=\"e2v\",\n", - " dimensions={\"lon\": \"glamv\", \"lat\": \"gphiv\"},\n", - " interp_method=\"nearest\",\n", - ")\n", - "fieldset.add_field(e1u)\n", - "fieldset.add_field(e2u)\n", - "fieldset.add_field(e1v)\n", - "fieldset.add_field(e2v)" + "field = fieldset.U.data[0, 0, :, :]\n", + "field = field.where(field != 0, np.nan) # Mask land values for better plotting\n", + "plt.pcolormesh(fieldset.U.grid.lon, fieldset.U.grid.lat, field, cmap=\"RdBu\")\n", + "\n", + "ds_out = xr.open_zarr(\"output_nemo3D.zarr\")\n", + "plt.scatter(ds_out.lon.T, ds_out.lat.T, c=ds_out.z.T, marker=\".\")\n", + "plt.show()" ] } ], diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index 1c91dd606..cca43eee1 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -30,12 +30,12 @@ TODO: Add links to Reference API throughout :titlesonly: examples/explanation_grids.md examples/tutorial_nemo_curvilinear.ipynb +examples/tutorial_nemo_3D.ipynb examples/tutorial_unitconverters.ipynb examples/tutorial_nestedgrids.ipynb ``` - diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 61394e7fc..601be0fbe 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -247,8 +247,10 @@ def from_copernicusmarine(ds: xr.Dataset): def from_nemo(ds: xr.Dataset): ds = ds.copy() - ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) + ds = _maybe_create_z_dim(ds) ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) + ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) + ds = _drop_unused_dimensions_and_coords(ds, _NEMO_DIMENSION_NAMES) ds = _maybe_rename_coords(ds, _NEMO_AXIS_VARNAMES) ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_NAMES) ds = _set_coords(ds, _NEMO_DIMENSION_NAMES) @@ -280,7 +282,7 @@ def from_nemo(ds: xr.Dataset): sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), ), - vertical_dimensions=(sgrid.DimDimPadding("depth_center", "depth", sgrid.Padding.LOW),), + vertical_dimensions=(sgrid.DimDimPadding("z_center", "z", sgrid.Padding.LOW),), ).to_attrs(), ) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") @@ -492,7 +494,7 @@ def _format_calendar_error_message(field: Field, reference_datetime: TimeLike) - "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], } -_NEMO_DIMENSION_NAMES = ["x", "y", "depth", "time", "lon", "lat"] +_NEMO_DIMENSION_NAMES = ["x", "y", "z", "time", "lon", "lat", "depth"] _NEMO_AXIS_VARNAMES = { "X": "lon", @@ -505,9 +507,24 @@ def _format_calendar_error_message(field: Field, reference_datetime: TimeLike) - "glamf": "lon", "gphif": "lat", "time_counter": "time", + "depthw": "depth", + "uo": "U", + "vo": "V", + "wo": "W", } +def _maybe_create_z_dim(ds): + if "depthw" in ds.dims: + lenZ = ds.sizes["depthw"] + for var in ds.data_vars: + for depthname in ["depthu", "depthv", "depthw"]: + if depthname in ds[var].dims: + ds[var] = ds[var].expand_dims(dim={"z": np.arange(lenZ)}, axis=1) + ds[var] = ds[var].isel({depthname: 0}, drop=True) + return ds + + def _maybe_rename_coords(ds, AXIS_VARNAMES): try: for axis, [coord] in ds.cf.axes.items(): @@ -531,6 +548,16 @@ def _assign_dims_as_coords(ds, DIMENSION_NAMES): return ds +def _drop_unused_dimensions_and_coords(ds, DIMENSION_NAMES): + for dim in ds.dims: + if dim not in DIMENSION_NAMES: + ds = ds.drop_dims(dim, errors="ignore") + for coord in ds.coords: + if coord not in DIMENSION_NAMES: + ds = ds.drop_vars(coord, errors="ignore") + return ds + + def _set_coords(ds, DIMENSION_NAMES): for varname in DIMENSION_NAMES: if varname in ds and varname not in ds.coords: From 6769a2c515b942f09032d27c8aadd7ad984010ca Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 16:53:31 +0100 Subject: [PATCH 05/38] Fixing negating vertical velocity code --- src/parcels/_core/fieldset.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 601be0fbe..237e881b2 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -224,7 +224,7 @@ def from_copernicusmarine(ds: xr.Dataset): ds = _maybe_rename_coords(ds, _COPERNICUS_MARINE_AXIS_VARNAMES) if "W" in ds.data_vars: # Negate W to convert from up positive to down positive (as that's the direction of positive z) - ds["W"] -= ds["W"] + ds["W"].data *= -1 if "grid" in ds.cf.cf_roles: raise ValueError( @@ -267,7 +267,7 @@ def from_nemo(ds: xr.Dataset): if "W" in ds.data_vars: # Negate W to convert from up positive to down positive (as that's the direction of positive z) - ds["W"] -= ds["W"] + ds["W"].data *= -1 if "grid" in ds.cf.cf_roles: raise ValueError( "Dataset already has a 'grid' variable (according to cf_roles). Didn't expect there to be grid metadata on copernicusmarine datasets - please open an issue with more information about your dataset." From d9f7a527bcec2a4b9514707bb6854a4ebf3bed86 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 17:06:02 +0100 Subject: [PATCH 06/38] Adding docstring for from_nemo ANd expading comments on loading in the dataset --- .../examples_v3/tutorial_nemo_3D.ipynb | 17 ++++++++++------ src/parcels/_core/fieldset.py | 20 +++++++++++++++++++ 2 files changed, 31 insertions(+), 6 deletions(-) diff --git a/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb b/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb index 805a293a8..1541e4612 100644 --- a/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb +++ b/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb @@ -55,6 +55,7 @@ "source": [ "data_folder = parcels.download_example_dataset(\"NemoNorthSeaORCA025-N006_data\")\n", "\n", + "# Open field datasets (with minimal loading to speed up)\n", "ds_fields = xr.open_mfdataset(\n", " data_folder.glob(\"ORCA*.nc\"),\n", " data_vars=\"minimal\",\n", @@ -62,19 +63,23 @@ " compat=\"override\",\n", ")\n", "\n", - "ds_coord = xr.open_dataset(data_folder / \"coordinates.nc\", decode_times=False).isel(\n", - " time=0, drop=True\n", - ")\n", - "ds_fields = xr.merge([ds_fields, ds_coord[[\"glamf\", \"gphif\"]]])\n", + "# Open coordinates dataset (with no time decoding because xarray otherwise complains)\n", + "ds_coords = xr.open_dataset(data_folder / \"coordinates.nc\", decode_times=False)\n", + "\n", + "# Remove time dimension from coordinates (otherwise xarray complains)\n", + "ds_coords = ds_coords.isel(time=0, drop=True)\n", + "\n", + "# Combine field and coordinate datasets\n", + "ds = xr.merge([ds_fields, ds_coords[[\"glamf\", \"gphif\"]]])\n", "\n", - "fieldset = parcels.FieldSet.from_nemo(ds_fields)" + "fieldset = parcels.FieldSet.from_nemo(ds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The code below is an example of how to then create a 3D simulation with particles, starting on a meridional line through the North Sea from the mouth of the river Rhine at 1m depth, and advecting them through the North Sea using the `AdvectionRK2_3D`" + "The code below is an example of how to then create a 3D simulation with particles, starting on a meridional line through the North Sea from the mouth of the river Rhine at 100m depth, and advecting them through the North Sea using the `AdvectionRK2_3D`" ] }, { diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 237e881b2..ee3eef912 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -246,6 +246,26 @@ def from_copernicusmarine(ds: xr.Dataset): return FieldSet.from_sgrid_conventions(ds, mesh="spherical") def from_nemo(ds: xr.Dataset): + """Create a FieldSet from a xarray.Dataset from NEMO netcdf files. + + Parameters + ---------- + ds : xarray.Dataset + xarray.Dataset as obtained from a set of NEMO netcdf files. + + Returns + ------- + FieldSet + FieldSet object containing the fields from the dataset that can be used for a Parcels simulation. + + Notes + ----- + The NEMO model (https://www.nemo-ocean.eu/) is used by a variety of oceanographic institutions around the world. + Output from these models may differ subtly in terms of variable names and metadata conventions. + This function attempts to standardize these differences to create a Parcels FieldSet. + If you encounter issues with your specific NEMO dataset, please open an issue on the Parcels GitHub repository with details about your dataset. + + """ ds = ds.copy() ds = _maybe_create_z_dim(ds) ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) From 1fa20b5c1fa75a146ad92275b9d06651ac831ad3 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 17:15:03 +0100 Subject: [PATCH 07/38] Support for NEMO fields without depth --- .../examples/tutorial_nemo_curvilinear.ipynb | 5 ----- src/parcels/_core/fieldset.py | 16 ++++++++++++++-- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index 6c980f03c..ca8c6aea3 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb @@ -74,11 +74,6 @@ " ds_fields.isel(time=0, drop=True).isel(z_a=0, drop=True).isel(z=0, drop=True)\n", ")\n", "\n", - "# TODO support 2D spatial fields in NEMO\n", - "ds_fields = ds_fields.expand_dims({\"depthw\": [0.0]})\n", - "ds_fields[\"glamf\"] = ds_fields[\"glamf\"].isel(depthw=0, drop=True)\n", - "ds_fields[\"gphif\"] = ds_fields[\"gphif\"].isel(depthw=0, drop=True)\n", - "\n", "fieldset = parcels.FieldSet.from_nemo(ds_fields)" ] }, diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index ee3eef912..768212ec3 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -267,13 +267,14 @@ def from_nemo(ds: xr.Dataset): """ ds = ds.copy() - ds = _maybe_create_z_dim(ds) + ds = _create_z_dim(ds) ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) ds = _drop_unused_dimensions_and_coords(ds, _NEMO_DIMENSION_NAMES) ds = _maybe_rename_coords(ds, _NEMO_AXIS_VARNAMES) ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_NAMES) ds = _set_coords(ds, _NEMO_DIMENSION_NAMES) + ds = _maybe_remove_z_from_lonlat(ds) ds = _set_axis_attrs(ds, _NEMO_AXIS_VARNAMES) expected_axes = set("XYZT") # TODO: Update after we have support for 2D spatial fields @@ -534,7 +535,7 @@ def _format_calendar_error_message(field: Field, reference_datetime: TimeLike) - } -def _maybe_create_z_dim(ds): +def _create_z_dim(ds): if "depthw" in ds.dims: lenZ = ds.sizes["depthw"] for var in ds.data_vars: @@ -542,6 +543,10 @@ def _maybe_create_z_dim(ds): if depthname in ds[var].dims: ds[var] = ds[var].expand_dims(dim={"z": np.arange(lenZ)}, axis=1) ds[var] = ds[var].isel({depthname: 0}, drop=True) + else: + if "z" not in ds.dims: + ds = ds.expand_dims({"z": [0]}) + ds["depth"] = xr.DataArray([0], dims=["z"]) return ds @@ -585,6 +590,13 @@ def _set_coords(ds, DIMENSION_NAMES): return ds +def _maybe_remove_z_from_lonlat(ds): + for coord in ["lon", "lat"]: + if coord in ds.coords and "z" in ds[coord].dims: + ds[coord] = ds[coord].squeeze("z", drop=True) + return ds + + def _set_axis_attrs(ds, AXIS_VARNAMES): for axis, varname in AXIS_VARNAMES.items(): if varname in ds.coords: From c0f5d3e1d2503043562287931cff19919d87da03 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 30 Dec 2025 17:20:03 +0100 Subject: [PATCH 08/38] Simplifying nemo dataset creation in tutorial --- docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index ca8c6aea3..3a174dd8e 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb @@ -70,9 +70,9 @@ " coords=\"minimal\",\n", " compat=\"override\",\n", ")\n", - "ds_fields = (\n", - " ds_fields.isel(time=0, drop=True).isel(z_a=0, drop=True).isel(z=0, drop=True)\n", - ")\n", + "\n", + "# Drop some auxiliary dimensions - otherwise Parcels will complain\n", + "ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True)\n", "\n", "fieldset = parcels.FieldSet.from_nemo(ds_fields)" ] From c5ad537aad23b0a14e34c7b33beb3918f9c3ef41 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 2 Jan 2026 13:45:09 +0100 Subject: [PATCH 09/38] Using glamf, gphif and depth as dimensions in from_nemo --- src/parcels/_core/fieldset.py | 56 +++++++++-------- tests/test_advection.py | 113 +++++----------------------------- 2 files changed, 48 insertions(+), 121 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 768212ec3..2892779a1 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -10,6 +10,7 @@ import xarray as xr import xgcm +from parcels._core.converters import GeographicPolar from parcels._core.field import Field, VectorField from parcels._core.utils import sgrid from parcels._core.utils.string import _assert_str_and_python_varname @@ -267,14 +268,15 @@ def from_nemo(ds: xr.Dataset): """ ds = ds.copy() - ds = _create_z_dim(ds) ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) + ds = _maybe_create_depth_dim(ds) + ds = _maybe_bring_UV_depths_to_depth(ds) ds = _drop_unused_dimensions_and_coords(ds, _NEMO_DIMENSION_NAMES) ds = _maybe_rename_coords(ds, _NEMO_AXIS_VARNAMES) ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_NAMES) ds = _set_coords(ds, _NEMO_DIMENSION_NAMES) - ds = _maybe_remove_z_from_lonlat(ds) + ds = _maybe_remove_depth_from_lonlat(ds) ds = _set_axis_attrs(ds, _NEMO_AXIS_VARNAMES) expected_axes = set("XYZT") # TODO: Update after we have support for 2D spatial fields @@ -298,15 +300,16 @@ def from_nemo(ds: xr.Dataset): attrs=sgrid.Grid2DMetadata( cf_role="grid_topology", topology_dimension=2, - node_dimensions=("lon", "lat"), + node_dimensions=("glamf", "gphif"), face_dimensions=( sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), ), - vertical_dimensions=(sgrid.DimDimPadding("z_center", "z", sgrid.Padding.LOW),), + vertical_dimensions=(sgrid.DimDimPadding("z_center", "depth", sgrid.Padding.HIGH),), ).to_attrs(), ) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") + fieldset.V.units = GeographicPolar() if "UV" in fieldset.fields: fieldset.UV.vector_interp_method = CGrid_Velocity if "UVW" in fieldset.fields: @@ -408,6 +411,11 @@ def from_sgrid_conventions( if "T" not in xgcm_kwargs["coords"]: xgcm_kwargs["coords"]["T"] = {"center": "time"} + if "lon" not in ds.coords or "lat" not in ds.coords: + node_dimensions = sgrid.load_mappings(ds.grid.node_dimensions) + ds["lon"] = ds[node_dimensions[0]] + ds["lat"] = ds[node_dimensions[1]] + # Create xgcm Grid object xgcm_grid = xgcm.Grid(ds, autoparse_metadata=False, **xgcm_kwargs, **_DEFAULT_XGCM_KWARGS) @@ -515,18 +523,16 @@ def _format_calendar_error_message(field: Field, reference_datetime: TimeLike) - "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], } -_NEMO_DIMENSION_NAMES = ["x", "y", "z", "time", "lon", "lat", "depth"] +_NEMO_DIMENSION_NAMES = ["x", "y", "time", "glamf", "gphif", "depth"] _NEMO_AXIS_VARNAMES = { - "X": "lon", - "Y": "lat", + "X": "glamf", + "Y": "gphif", "Z": "depth", "T": "time", } _NEMO_VARNAMES_MAPPING = { - "glamf": "lon", - "gphif": "lat", "time_counter": "time", "depthw": "depth", "uo": "U", @@ -535,18 +541,18 @@ def _format_calendar_error_message(field: Field, reference_datetime: TimeLike) - } -def _create_z_dim(ds): - if "depthw" in ds.dims: - lenZ = ds.sizes["depthw"] - for var in ds.data_vars: - for depthname in ["depthu", "depthv", "depthw"]: - if depthname in ds[var].dims: - ds[var] = ds[var].expand_dims(dim={"z": np.arange(lenZ)}, axis=1) - ds[var] = ds[var].isel({depthname: 0}, drop=True) - else: - if "z" not in ds.dims: - ds = ds.expand_dims({"z": [0]}) - ds["depth"] = xr.DataArray([0], dims=["z"]) +def _maybe_bring_UV_depths_to_depth(ds): + if "U" in ds.variables and "depthu" in ds.U.coords and "depth" in ds.coords: + ds["U"] = ds["U"].assign_coords(depthu=ds["depth"].values).rename({"depthu": "depth"}) + if "V" in ds.variables and "depthv" in ds.V.coords and "depth" in ds.coords: + ds["V"] = ds["V"].assign_coords(depthv=ds["depth"].values).rename({"depthv": "depth"}) + return ds + + +def _maybe_create_depth_dim(ds): + if "depth" not in ds.dims: + ds = ds.expand_dims({"depth": [0]}) + ds["depth"] = xr.DataArray([0], dims=["depth"]) return ds @@ -590,10 +596,10 @@ def _set_coords(ds, DIMENSION_NAMES): return ds -def _maybe_remove_z_from_lonlat(ds): - for coord in ["lon", "lat"]: - if coord in ds.coords and "z" in ds[coord].dims: - ds[coord] = ds[coord].squeeze("z", drop=True) +def _maybe_remove_depth_from_lonlat(ds): + for coord in ["glamf", "gphif"]: + if coord in ds.coords and "depth" in ds[coord].dims: + ds[coord] = ds[coord].squeeze("depth", drop=True) return ds diff --git a/tests/test_advection.py b/tests/test_advection.py index 0a0b7484a..2a90cec4f 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -1,7 +1,6 @@ import numpy as np import pytest import xarray as xr -import xgcm import parcels from parcels import Field, FieldSet, Particle, ParticleFile, ParticleSet, StatusCode, Variable, VectorField, XGrid @@ -417,32 +416,14 @@ def UpdateP(particles, fieldset): # pragma: no cover def test_nemo_curvilinear_fieldset(): data_folder = parcels.download_example_dataset("NemoCurvilinear_data") - files = data_folder.glob("*.nc4") - ds = xr.open_mfdataset(files, combine="nested", data_vars="minimal", coords="minimal", compat="override") - ds = ( - ds.isel(time_counter=0, drop=True) - .isel(time=0, drop=True) - .isel(z_a=0, drop=True) - .rename({"glamf": "lon", "gphif": "lat", "z": "depth"}) + ds_fields = xr.open_mfdataset( + data_folder.glob("*.nc4"), + data_vars="minimal", + coords="minimal", + compat="override", ) - - xgcm_grid = xgcm.Grid( - ds, - coords={ - "X": {"left": "x"}, - "Y": {"left": "y"}, - }, - periodic=False, - autoparse_metadata=False, - ) - grid = XGrid(xgcm_grid, mesh="spherical") - - U = parcels.Field("U", ds["U"], grid, interp_method=XLinear) - V = parcels.Field("V", ds["V"], grid, interp_method=XLinear) - U.units = parcels.GeographicPolar() - V.units = parcels.GeographicPolar() # U and V need GeographicPolar for C-Grid interpolation to work correctly - UV = parcels.VectorField("UV", U, V, vector_interp_method=CGrid_Velocity) - fieldset = parcels.FieldSet([U, V, UV]) + ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True) + fieldset = parcels.FieldSet.from_nemo(ds_fields) npart = 20 lonp = 30 * np.ones(npart) @@ -456,77 +437,17 @@ def test_nemo_curvilinear_fieldset(): @pytest.mark.parametrize("kernel", [AdvectionRK4, AdvectionRK4_3D]) def test_nemo_3D_curvilinear_fieldset(kernel): - download_dir = parcels.download_example_dataset("NemoNorthSeaORCA025-N006_data") - ufiles = download_dir.glob("*U.nc") - dsu = xr.open_mfdataset(ufiles, decode_times=False, drop_variables=["nav_lat", "nav_lon"]) - dsu = dsu.rename({"time_counter": "time", "uo": "U"}) - - vfiles = download_dir.glob("*V.nc") - dsv = xr.open_mfdataset(vfiles, decode_times=False, drop_variables=["nav_lat", "nav_lon"]) - dsv = dsv.rename({"time_counter": "time", "vo": "V"}) - - wfiles = download_dir.glob("*W.nc") - dsw = xr.open_mfdataset(wfiles, decode_times=False, drop_variables=["nav_lat", "nav_lon"]) - dsw = dsw.rename({"time_counter": "time", "depthw": "depth", "wo": "W"}) - - dsu = dsu.assign_coords(depthu=dsw.depth.values) - dsu = dsu.rename({"depthu": "depth"}) - - dsv = dsv.assign_coords(depthv=dsw.depth.values) - dsv = dsv.rename({"depthv": "depth"}) - - coord_file = f"{download_dir}/coordinates.nc" - dscoord = xr.open_dataset(coord_file, decode_times=False).rename({"glamf": "lon", "gphif": "lat"}) - dscoord = dscoord.isel(time=0, drop=True) - - ds = xr.merge([dsu, dsv, dsw, dscoord]) - ds = ds.drop_vars( - [ - "uos", - "vos", - "nav_lev", - "nav_lon", - "nav_lat", - "tauvo", - "tauuo", - "time_steps", - "gphiu", - "gphiv", - "gphit", - "glamu", - "glamv", - "glamt", - "time_centered_bounds", - "time_counter_bounds", - "time_centered", - ] + data_folder = parcels.download_example_dataset("NemoNorthSeaORCA025-N006_data") + ds_fields = xr.open_mfdataset( + data_folder.glob("ORCA*.nc"), + data_vars="minimal", + coords="minimal", + compat="override", ) - ds = ds.drop_vars(["e1f", "e1t", "e1u", "e1v", "e2f", "e2t", "e2u", "e2v"]) - ds["time"] = [np.timedelta64(int(t), "s") + np.datetime64("1900-01-01") for t in ds["time"]] - - ds["W"] *= -1 # Invert W velocity - - xgcm_grid = xgcm.Grid( - ds, - coords={ - "X": {"left": "x"}, - "Y": {"left": "y"}, - "Z": {"left": "depth"}, - "T": {"center": "time"}, - }, - periodic=False, - autoparse_metadata=False, - ) - grid = XGrid(xgcm_grid, mesh="spherical") - - U = parcels.Field("U", ds["U"], grid, interp_method=XLinear) - V = parcels.Field("V", ds["V"], grid, interp_method=XLinear) - W = parcels.Field("W", ds["W"], grid, interp_method=XLinear) - U.units = parcels.GeographicPolar() - V.units = parcels.GeographicPolar() # U and V need GoegraphicPolar for C-Grid interpolation to work correctly - UV = parcels.VectorField("UV", U, V, vector_interp_method=CGrid_Velocity) - UVW = parcels.VectorField("UVW", U, V, W, vector_interp_method=CGrid_Velocity) - fieldset = parcels.FieldSet([U, V, W, UV, UVW]) + ds_coords = xr.open_dataset(data_folder / "coordinates.nc", decode_times=False) + ds_coords = ds_coords.isel(time=0, drop=True) + ds = xr.merge([ds_fields, ds_coords[["glamf", "gphif"]]]) + fieldset = parcels.FieldSet.from_nemo(ds) npart = 10 lons = np.linspace(1.9, 3.4, npart) From 60e2f387bfa8be13d869c49a92080ba2b89bca00 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 2 Jan 2026 17:18:51 +0100 Subject: [PATCH 10/38] Using xgcm coordinates for cgrid interpolation selection --- src/parcels/interpolators.py | 64 ++++++++++++++++-------------------- 1 file changed, 28 insertions(+), 36 deletions(-) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index 8e23d4952..3d0dc91df 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -237,33 +237,20 @@ def CGrid_Velocity( corner_data = corner_data[0, :, :] * (1 - tau_full) + corner_data[1, :, :] * tau_full else: corner_data = corner_data[0, :, :] - # # See code below for v3 version - # # if self.gridindexingtype == "nemo": - # # U0 = self.U.data[ti, zi, yi + 1, xi] * c4 - # # U1 = self.U.data[ti, zi, yi + 1, xi + 1] * c2 - # # V0 = self.V.data[ti, zi, yi, xi + 1] * c1 - # # V1 = self.V.data[ti, zi, yi + 1, xi + 1] * c3 - # # elif self.gridindexingtype in ["mitgcm", "croco"]: - # # U0 = self.U.data[ti, zi, yi, xi] * c4 - # # U1 = self.U.data[ti, zi, yi, xi + 1] * c2 - # # V0 = self.V.data[ti, zi, yi, xi] * c1 - # # V1 = self.V.data[ti, zi, yi + 1, xi] * c3 - # # TODO Nick can you help use xgcm to fix this implementation? - - # # CROCO and MITgcm grid indexing, - # if data is U: - # U0 = corner_data[:, 0] * c4 - # U1 = corner_data[:, 1] * c2 - # elif data is V: - # V0 = corner_data[:, 0] * c1 - # V1 = corner_data[:, 2] * c3 - # # NEMO grid indexing + + X_axis_coord = vectorfield.U.grid.xgcm_grid.axes["X"].coords.keys() + Y_axis_coord = vectorfield.V.grid.xgcm_grid.axes["Y"].coords.keys() + if data is U: - U0 = corner_data[:, 2] * c4 - U1 = corner_data[:, 3] * c2 + val0 = 2 if "right" in X_axis_coord else 0 + val1 = 3 if "right" in X_axis_coord else 1 + U0 = corner_data[:, val0] * c4 + U1 = corner_data[:, val1] * c2 elif data is V: - V0 = corner_data[:, 1] * c1 - V1 = corner_data[:, 3] * c3 + val0 = 1 if "right" in Y_axis_coord else 0 + val1 = 3 if "right" in Y_axis_coord else 2 + V0 = corner_data[:, val0] * c1 + V1 = corner_data[:, val1] * c3 U = (1 - xsi) * U0 + xsi * U1 V = (1 - eta) * V0 + eta * V1 @@ -304,17 +291,22 @@ def CGrid_Velocity( ti_1 = np.clip(ti + 1, 0, tdim - 1) ti_full = np.concatenate([np.repeat(ti, 2), np.repeat(ti_1, 2)]) - # Z coordinates: 1 points at zi, repeated for both time levels - zi_1 = np.clip(zi + 1, 0, zdim - 1) - zi_full = np.tile(np.array([zi, zi_1]).flatten(), lenT) - - # Y coordinates: yi+1 for each spatial point, repeated for time/z - yi_1 = np.clip(yi + 1, 0, ydim - 1) - yi_full = np.tile(yi_1, (lenT) * 2) - - # X coordinates: xi+1 for each spatial point, repeated for time/z - xi_1 = np.clip(xi + 1, 0, xdim - 1) - xi_full = np.tile(xi_1, (lenT) * 2) + # Z coordinates: 1 points at zi + offset, repeated for both time levels + Z_axis_coord = vectorfield.W.grid.xgcm_grid.axes["Z"].coords.keys() + offset = 1 if "right" in Z_axis_coord else 0 + zi_0 = np.clip(zi + offset, 0, zdim - 1) + zi_1 = np.clip(zi + 1 + offset, 0, zdim - 1) + zi_full = np.tile(np.array([zi_0, zi_1]).flatten(), lenT) + + # Y coordinates: yi+offset for each spatial point, repeated for time/z + offset = 1 if "right" in Y_axis_coord else 0 + yi_0 = np.clip(yi + offset, 0, ydim - 1) + yi_full = np.tile(yi_0, (lenT) * 2) + + # X coordinates: xi+offset for each spatial point, repeated for time/z + offset = 1 if "right" in X_axis_coord else 0 + xi_o = np.clip(xi + offset, 0, xdim - 1) + xi_full = np.tile(xi_o, (lenT) * 2) axis_dim = grid.get_axis_dim_mapping(data.dims) From d36cf8270b08e17f31a743866241e490b6879bfd Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 5 Jan 2026 10:42:50 +0100 Subject: [PATCH 11/38] Reordering corner_data in CGrid_Velocity To use same ordering as in XLinear --- src/parcels/interpolators.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index 3d0dc91df..3e93ccb8c 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -208,14 +208,14 @@ def CGrid_Velocity( # Y coordinates: [yi, yi, yi+1, yi+1] for each spatial point, repeated for time/z yi_1 = np.clip(yi + 1, 0, ydim - 1) - yi_full = np.tile(np.repeat(np.column_stack([yi, yi_1]), 2), (lenT)) + yi_full = np.tile(np.array([yi, yi, yi_1, yi_1]).flatten(), lenT) # # TODO check why in some cases minus needed here!!! # yi_minus_1 = np.clip(yi - 1, 0, ydim - 1) # yi = np.tile(np.repeat(np.column_stack([yi_minus_1, yi]), 2), (lenT)) # X coordinates: [xi, xi+1, xi, xi+1] for each spatial point, repeated for time/z xi_1 = np.clip(xi + 1, 0, xdim - 1) - xi_full = np.tile(np.column_stack([xi, xi_1, xi, xi_1]).flatten(), (lenT)) + xi_full = np.tile(np.array([xi, xi_1]).flatten(), lenT * 2) for data in [U, V]: axis_dim = grid.get_axis_dim_mapping(data.dims) @@ -230,27 +230,25 @@ def CGrid_Velocity( if "time" in data.dims: selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) - corner_data = data.isel(selection_dict).data.reshape(lenT, len(xsi), 4) + corner_data = data.isel(selection_dict).data.reshape(lenT, 2, 2, len(xsi)) if lenT == 2: - tau_full = tau[:, np.newaxis] - corner_data = corner_data[0, :, :] * (1 - tau_full) + corner_data[1, :, :] * tau_full + tau_full = tau[np.newaxis, :] + corner_data = corner_data[0, :] * (1 - tau_full) + corner_data[1, :] * tau_full else: - corner_data = corner_data[0, :, :] + corner_data = corner_data[0, :] X_axis_coord = vectorfield.U.grid.xgcm_grid.axes["X"].coords.keys() Y_axis_coord = vectorfield.V.grid.xgcm_grid.axes["Y"].coords.keys() if data is U: - val0 = 2 if "right" in X_axis_coord else 0 - val1 = 3 if "right" in X_axis_coord else 1 - U0 = corner_data[:, val0] * c4 - U1 = corner_data[:, val1] * c2 + offset = 1 if "right" in X_axis_coord else 0 + U0 = corner_data[offset, 0, :] * c4 + U1 = corner_data[offset, 1, :] * c2 elif data is V: - val0 = 1 if "right" in Y_axis_coord else 0 - val1 = 3 if "right" in Y_axis_coord else 2 - V0 = corner_data[:, val0] * c1 - V1 = corner_data[:, val1] * c3 + offset = 1 if "right" in Y_axis_coord else 0 + V0 = corner_data[0, offset, :] * c1 + V1 = corner_data[1, offset, :] * c3 U = (1 - xsi) * U0 + xsi * U1 V = (1 - eta) * V0 + eta * V1 From 0f91c235c3a1706676cc295bce942a6ec2452a94 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 5 Jan 2026 12:08:23 +0100 Subject: [PATCH 12/38] Separating U and V interpolation in CGrid_Velocity This requires smaller selection_dict for the isel, so hopefully faster code --- src/parcels/interpolators.py | 124 ++++++++++++++++++++--------------- 1 file changed, 71 insertions(+), 53 deletions(-) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index 3e93ccb8c..ef115e63a 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -193,81 +193,99 @@ def CGrid_Velocity( lenT = 2 if np.any(tau > 0) else 1 - # Create arrays of corner points for xarray.isel - # TODO C grid may not need all xi and yi cornerpoints, so could speed up here? - - # Time coordinates: 4 points at ti, then 4 points at ti+1 + # Time coordinates: 2 points at ti, then 2 points at ti+1 if lenT == 1: - ti_full = np.repeat(ti, 4) + ti_full = np.repeat(ti, 2) else: ti_1 = np.clip(ti + 1, 0, tdim - 1) - ti_full = np.concatenate([np.repeat(ti, 4), np.repeat(ti_1, 4)]) + ti_full = np.concatenate([np.repeat(ti, 2), np.repeat(ti_1, 2)]) - # Z coordinates: 4 points at zi, repeated for both time levels - zi_full = np.repeat(zi, lenT * 4) + # Z coordinates: 2 points at zi, repeated for both time levels + zi_full = np.repeat(zi, lenT * 2) - # Y coordinates: [yi, yi, yi+1, yi+1] for each spatial point, repeated for time/z - yi_1 = np.clip(yi + 1, 0, ydim - 1) - yi_full = np.tile(np.array([yi, yi, yi_1, yi_1]).flatten(), lenT) - # # TODO check why in some cases minus needed here!!! - # yi_minus_1 = np.clip(yi - 1, 0, ydim - 1) - # yi = np.tile(np.repeat(np.column_stack([yi_minus_1, yi]), 2), (lenT)) + # Y coordinates for U: [yi+offset, yi+offset] for each spatial point, repeated for time/z + Y_axis_coord = vectorfield.U.grid.xgcm_grid.axes["Y"].coords.keys() + offset = 1 if "right" in Y_axis_coord else 0 + yi_o = np.clip(yi + offset, 0, ydim - 1) + yi_full = np.tile(np.array([yi_o, yi_o]).flatten(), lenT) - # X coordinates: [xi, xi+1, xi, xi+1] for each spatial point, repeated for time/z + # X coordinates for U: [xi, xi+1] for each spatial point, repeated for time/z xi_1 = np.clip(xi + 1, 0, xdim - 1) - xi_full = np.tile(np.array([xi, xi_1]).flatten(), lenT * 2) + xi_full = np.tile(np.array([xi, xi_1]).flatten(), lenT) - for data in [U, V]: - axis_dim = grid.get_axis_dim_mapping(data.dims) + axis_dim = grid.get_axis_dim_mapping(U.dims) - # Create DataArrays for indexing - selection_dict = { - axis_dim["X"]: xr.DataArray(xi_full, dims=("points")), - axis_dim["Y"]: xr.DataArray(yi_full, dims=("points")), - } - if "Z" in axis_dim: - selection_dict[axis_dim["Z"]] = xr.DataArray(zi_full, dims=("points")) - if "time" in data.dims: - selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) + # Create DataArrays for indexing + selection_dict = { + axis_dim["X"]: xr.DataArray(xi_full, dims=("points")), + axis_dim["Y"]: xr.DataArray(yi_full, dims=("points")), + } + if "Z" in axis_dim: + selection_dict[axis_dim["Z"]] = xr.DataArray(zi_full, dims=("points")) + if "time" in U.dims: + selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) - corner_data = data.isel(selection_dict).data.reshape(lenT, 2, 2, len(xsi)) + corner_data = U.isel(selection_dict).data.reshape(lenT, 2, len(xsi)) - if lenT == 2: - tau_full = tau[np.newaxis, :] - corner_data = corner_data[0, :] * (1 - tau_full) + corner_data[1, :] * tau_full - else: - corner_data = corner_data[0, :] + if lenT == 2: + tau_full = tau[np.newaxis, :] + corner_data = corner_data[0, :] * (1 - tau_full) + corner_data[1, :] * tau_full + else: + corner_data = corner_data[0, :] + + U0 = corner_data[0, :] * c4 + U1 = corner_data[1, :] * c2 + Uvel = (1 - xsi) * U0 + xsi * U1 - X_axis_coord = vectorfield.U.grid.xgcm_grid.axes["X"].coords.keys() - Y_axis_coord = vectorfield.V.grid.xgcm_grid.axes["Y"].coords.keys() + # Y coordinates for V: [yi, yi+1] for each spatial point, repeated for time/z + yi_1 = np.clip(yi + 1, 0, ydim - 1) + yi_full = np.tile(np.array([yi, yi_1]).flatten(), lenT) - if data is U: - offset = 1 if "right" in X_axis_coord else 0 - U0 = corner_data[offset, 0, :] * c4 - U1 = corner_data[offset, 1, :] * c2 - elif data is V: - offset = 1 if "right" in Y_axis_coord else 0 - V0 = corner_data[0, offset, :] * c1 - V1 = corner_data[1, offset, :] * c3 + # X coordinates for V: [xi+offset, xi+offset] for each spatial point, repeated for time/z + X_axis_coord = vectorfield.V.grid.xgcm_grid.axes["X"].coords.keys() + offset = 1 if "right" in X_axis_coord else 0 + xi_o = np.clip(xi + offset, 0, xdim - 1) + xi_full = np.tile(np.array([xi_o, xi_o]).flatten(), lenT) + + axis_dim = grid.get_axis_dim_mapping(V.dims) + + # Create DataArrays for indexing + selection_dict = { + axis_dim["X"]: xr.DataArray(xi_full, dims=("points")), + axis_dim["Y"]: xr.DataArray(yi_full, dims=("points")), + } + if "Z" in axis_dim: + selection_dict[axis_dim["Z"]] = xr.DataArray(zi_full, dims=("points")) + if "time" in V.dims: + selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) + + corner_data = V.isel(selection_dict).data.reshape(lenT, 2, len(xsi)) + + if lenT == 2: + tau_full = tau[np.newaxis, :] + corner_data = corner_data[0, :] * (1 - tau_full) + corner_data[1, :] * tau_full + else: + corner_data = corner_data[0, :] - U = (1 - xsi) * U0 + xsi * U1 - V = (1 - eta) * V0 + eta * V1 + V0 = corner_data[0, :] * c1 + V1 = corner_data[1, :] * c3 + Vvel = (1 - eta) * V0 + eta * V1 meshJac = 1852 * 60.0 if grid._mesh == "spherical" else 1 jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac u = ( - (-(1 - eta) * U - (1 - xsi) * V) * px[0] - + ((1 - eta) * U - xsi * V) * px[1] - + (eta * U + xsi * V) * px[2] - + (-eta * U + (1 - xsi) * V) * px[3] + (-(1 - eta) * Uvel - (1 - xsi) * Vvel) * px[0] + + ((1 - eta) * Uvel - xsi * Vvel) * px[1] + + (eta * Uvel + xsi * Vvel) * px[2] + + (-eta * Uvel + (1 - xsi) * Vvel) * px[3] ) / jac v = ( - (-(1 - eta) * U - (1 - xsi) * V) * py[0] - + ((1 - eta) * U - xsi * V) * py[1] - + (eta * U + xsi * V) * py[2] - + (-eta * U + (1 - xsi) * V) * py[3] + (-(1 - eta) * Uvel - (1 - xsi) * Vvel) * py[0] + + ((1 - eta) * Uvel - xsi * Vvel) * py[1] + + (eta * Uvel + xsi * Vvel) * py[2] + + (-eta * Uvel + (1 - xsi) * Vvel) * py[3] ) / jac if is_dask_collection(u): u = u.compute() From 17c9d4e2f17d292d2afe3b44a5fd112763d86b4a Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 5 Jan 2026 13:30:17 +0100 Subject: [PATCH 13/38] Cleaning up CGrid_Velocity code --- .../examples_v3/tutorial_nemo_3D.ipynb | 3 +- src/parcels/interpolators.py | 162 +++++++----------- 2 files changed, 65 insertions(+), 100 deletions(-) diff --git a/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb b/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb index 1541e4612..a70f5a58e 100644 --- a/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb +++ b/docs/user_guide/examples_v3/tutorial_nemo_3D.ipynb @@ -130,7 +130,8 @@ "plt.pcolormesh(fieldset.U.grid.lon, fieldset.U.grid.lat, field, cmap=\"RdBu\")\n", "\n", "ds_out = xr.open_zarr(\"output_nemo3D.zarr\")\n", - "plt.scatter(ds_out.lon.T, ds_out.lat.T, c=ds_out.z.T, marker=\".\")\n", + "plt.scatter(ds_out.lon.T, ds_out.lat.T, c=-ds_out.z.T, marker=\".\")\n", + "plt.colorbar(label=\"Depth (m)\")\n", "plt.show()" ] } diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index ef115e63a..377f0d79a 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -166,6 +166,12 @@ def CGrid_Velocity( V = vectorfield.V.data grid = vectorfield.grid tdim, zdim, ydim, xdim = U.shape[0], U.shape[1], U.shape[2], U.shape[3] + lenT = 2 if np.any(tau > 0) else 1 + + offsets = {} + for axis in ["X", "Y"]: + axis_coords = grid.xgcm_grid.axes[axis].coords.keys() + offsets[axis] = 1 if "right" in axis_coords else 0 if grid.lon.ndim == 1: px = np.array([grid.lon[xi], grid.lon[xi + 1], grid.lon[xi + 1], grid.lon[xi]]) @@ -191,88 +197,77 @@ def CGrid_Velocity( py[3], py[0], px[3], px[0], grid._mesh, np.einsum("ij,ji->i", i_u.phi2D_lin(eta, 0.0), py) ) - lenT = 2 if np.any(tau > 0) else 1 + def _create_selection_dict(dims, zdir=False): + """Helper function to create DataArrays for indexing.""" + axis_dim = grid.get_axis_dim_mapping(dims) + selection_dict = { + axis_dim["X"]: xr.DataArray(xi_full, dims=("points")), + axis_dim["Y"]: xr.DataArray(yi_full, dims=("points")), + } - # Time coordinates: 2 points at ti, then 2 points at ti+1 - if lenT == 1: - ti_full = np.repeat(ti, 2) - else: - ti_1 = np.clip(ti + 1, 0, tdim - 1) - ti_full = np.concatenate([np.repeat(ti, 2), np.repeat(ti_1, 2)]) + # Time coordinates: 2 points at ti, then 2 points at ti+1 + if "time" in dims: + if lenT == 1: + ti_full = np.repeat(ti, 2) + else: + ti_1 = np.clip(ti + 1, 0, tdim - 1) + ti_full = np.concatenate([np.repeat(ti, 2), np.repeat(ti_1, 2)]) + selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) - # Z coordinates: 2 points at zi, repeated for both time levels - zi_full = np.repeat(zi, lenT * 2) + if "Z" in axis_dim: + if zdir: + # Z coordinates: 1 point at zi and 1 point at zi+1 repeated for lenT time levels + zi_0 = np.clip(zi, 0, zdim - 1) + zi_1 = np.clip(zi + 1, 0, zdim - 1) + zi_full = np.tile(np.array([zi_0, zi_1]).flatten(), lenT) + else: + # Z coordinates: 2 points at zi, repeated for lenT time levels + zi_full = np.repeat(zi, lenT * 2) + selection_dict[axis_dim["Z"]] = xr.DataArray(zi_full, dims=("points")) + + return selection_dict + + def _compute_corner_data(data, selection_dict) -> np.ndarray: + """Helper function to load and reduce corner data over time dimension if needed.""" + corner_data = data.isel(selection_dict).data.reshape(lenT, 2, len(xsi)) + + if lenT == 2: + tau_full = tau[np.newaxis, :] + corner_data = corner_data[0, :] * (1 - tau_full) + corner_data[1, :] * tau_full + else: + corner_data = corner_data[0, :] + return corner_data - # Y coordinates for U: [yi+offset, yi+offset] for each spatial point, repeated for time/z - Y_axis_coord = vectorfield.U.grid.xgcm_grid.axes["Y"].coords.keys() - offset = 1 if "right" in Y_axis_coord else 0 - yi_o = np.clip(yi + offset, 0, ydim - 1) + # Compute U velocity + yi_o = np.clip(yi + offsets["Y"], 0, ydim - 1) yi_full = np.tile(np.array([yi_o, yi_o]).flatten(), lenT) - # X coordinates for U: [xi, xi+1] for each spatial point, repeated for time/z xi_1 = np.clip(xi + 1, 0, xdim - 1) xi_full = np.tile(np.array([xi, xi_1]).flatten(), lenT) - axis_dim = grid.get_axis_dim_mapping(U.dims) - - # Create DataArrays for indexing - selection_dict = { - axis_dim["X"]: xr.DataArray(xi_full, dims=("points")), - axis_dim["Y"]: xr.DataArray(yi_full, dims=("points")), - } - if "Z" in axis_dim: - selection_dict[axis_dim["Z"]] = xr.DataArray(zi_full, dims=("points")) - if "time" in U.dims: - selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) - - corner_data = U.isel(selection_dict).data.reshape(lenT, 2, len(xsi)) - - if lenT == 2: - tau_full = tau[np.newaxis, :] - corner_data = corner_data[0, :] * (1 - tau_full) + corner_data[1, :] * tau_full - else: - corner_data = corner_data[0, :] + selection_dict = _create_selection_dict(U.dims) + corner_data = _compute_corner_data(U, selection_dict) U0 = corner_data[0, :] * c4 U1 = corner_data[1, :] * c2 Uvel = (1 - xsi) * U0 + xsi * U1 - # Y coordinates for V: [yi, yi+1] for each spatial point, repeated for time/z + # Compute V velocity yi_1 = np.clip(yi + 1, 0, ydim - 1) yi_full = np.tile(np.array([yi, yi_1]).flatten(), lenT) - # X coordinates for V: [xi+offset, xi+offset] for each spatial point, repeated for time/z - X_axis_coord = vectorfield.V.grid.xgcm_grid.axes["X"].coords.keys() - offset = 1 if "right" in X_axis_coord else 0 - xi_o = np.clip(xi + offset, 0, xdim - 1) + xi_o = np.clip(xi + offsets["X"], 0, xdim - 1) xi_full = np.tile(np.array([xi_o, xi_o]).flatten(), lenT) - axis_dim = grid.get_axis_dim_mapping(V.dims) - - # Create DataArrays for indexing - selection_dict = { - axis_dim["X"]: xr.DataArray(xi_full, dims=("points")), - axis_dim["Y"]: xr.DataArray(yi_full, dims=("points")), - } - if "Z" in axis_dim: - selection_dict[axis_dim["Z"]] = xr.DataArray(zi_full, dims=("points")) - if "time" in V.dims: - selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) - - corner_data = V.isel(selection_dict).data.reshape(lenT, 2, len(xsi)) - - if lenT == 2: - tau_full = tau[np.newaxis, :] - corner_data = corner_data[0, :] * (1 - tau_full) + corner_data[1, :] * tau_full - else: - corner_data = corner_data[0, :] + selection_dict = _create_selection_dict(V.dims) + corner_data = _compute_corner_data(V, selection_dict) V0 = corner_data[0, :] * c1 V1 = corner_data[1, :] * c3 Vvel = (1 - eta) * V0 + eta * V1 - meshJac = 1852 * 60.0 if grid._mesh == "spherical" else 1 - + # rotate velocities to eastward/northward directions + meshJac = 1852 * 60.0 if grid._mesh == "spherical" else 1.0 jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac u = ( @@ -299,49 +294,18 @@ def CGrid_Velocity( u = np.where(np.abs(dlon / lon) > 1e-4, np.nan, u) if vectorfield.W: - data = vectorfield.W.data - # Time coordinates: 2 points at ti, then 2 points at ti+1 - if lenT == 1: - ti_full = np.repeat(ti, 2) - else: - ti_1 = np.clip(ti + 1, 0, tdim - 1) - ti_full = np.concatenate([np.repeat(ti, 2), np.repeat(ti_1, 2)]) - - # Z coordinates: 1 points at zi + offset, repeated for both time levels - Z_axis_coord = vectorfield.W.grid.xgcm_grid.axes["Z"].coords.keys() - offset = 1 if "right" in Z_axis_coord else 0 - zi_0 = np.clip(zi + offset, 0, zdim - 1) - zi_1 = np.clip(zi + 1 + offset, 0, zdim - 1) - zi_full = np.tile(np.array([zi_0, zi_1]).flatten(), lenT) - - # Y coordinates: yi+offset for each spatial point, repeated for time/z - offset = 1 if "right" in Y_axis_coord else 0 - yi_0 = np.clip(yi + offset, 0, ydim - 1) - yi_full = np.tile(yi_0, (lenT) * 2) - - # X coordinates: xi+offset for each spatial point, repeated for time/z - offset = 1 if "right" in X_axis_coord else 0 - xi_o = np.clip(xi + offset, 0, xdim - 1) - xi_full = np.tile(xi_o, (lenT) * 2) - - axis_dim = grid.get_axis_dim_mapping(data.dims) + W = vectorfield.W.data - # Create DataArrays for indexing - selection_dict = { - axis_dim["X"]: xr.DataArray(xi_full, dims=("points")), - axis_dim["Y"]: xr.DataArray(yi_full, dims=("points")), - axis_dim["Z"]: xr.DataArray(zi_full, dims=("points")), - } - if "time" in data.dims: - selection_dict["time"] = xr.DataArray(ti_full, dims=("points")) + # Y coordinates: yi+offset for each spatial point, repeated for time + yi_o = np.clip(yi + offsets["Y"], 0, ydim - 1) + yi_full = np.tile(yi_o, (lenT) * 2) - corner_data = data.isel(selection_dict).data.reshape(lenT, 2, len(xsi)) + # X coordinates: xi+offset for each spatial point, repeated for time + xi_o = np.clip(xi + offsets["X"], 0, xdim - 1) + xi_full = np.tile(xi_o, (lenT) * 2) - if lenT == 2: - tau_full = tau[np.newaxis, :] - corner_data = corner_data[0, :, :] * (1 - tau_full) + corner_data[1, :, :] * tau_full - else: - corner_data = corner_data[0, :, :] + selection_dict = _create_selection_dict(W.dims, zdir=True) + corner_data = _compute_corner_data(W, selection_dict) w = corner_data[0, :] * (1 - zeta) + corner_data[1, :] * zeta if is_dask_collection(w): From 79b4c3788b3ccdf278c0fe9367d843970b6ddf79 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 6 Jan 2026 11:38:28 +0100 Subject: [PATCH 14/38] Fixing bug in CGrid_Tracer interpolation --- src/parcels/interpolators.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index 377f0d79a..f40c6d90e 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -340,12 +340,9 @@ def CGrid_Tracer( if lenT == 2: ti_1 = np.clip(ti + 1, 0, data.shape[0] - 1) ti = np.concatenate([np.repeat(ti), np.repeat(ti_1)]) - zi_1 = np.clip(zi + 1, 0, data.shape[1] - 1) - zi = np.concatenate([np.repeat(zi), np.repeat(zi_1)]) - yi_1 = np.clip(yi + 1, 0, data.shape[2] - 1) - yi = np.concatenate([np.repeat(yi), np.repeat(yi_1)]) - xi_1 = np.clip(xi + 1, 0, data.shape[3] - 1) - xi = np.concatenate([np.repeat(xi), np.repeat(xi_1)]) + zi = np.tile(zi, (lenT) * 2) + yi = np.tile(yi, (lenT) * 2) + xi = np.tile(xi, (lenT) * 2) # Create DataArrays for indexing selection_dict = { From 7343a437344efdca0469e246e588e0686a3fb669 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 6 Jan 2026 11:49:40 +0100 Subject: [PATCH 15/38] Using offsets in CGrid_Tracer And also separating offset calculation into its own helper function --- src/parcels/interpolators.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index f40c6d90e..87412d00d 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -97,6 +97,19 @@ def _get_corner_data_Agrid( return data.isel(selection_dict).data.reshape(lenT, lenZ, 2, 2, npart) +def _get_offsets_dictionary(grid): + offsets = {} + for axis in ["X", "Y"]: + axis_coords = grid.xgcm_grid.axes[axis].coords.keys() + offsets[axis] = 1 if "right" in axis_coords else 0 + if "Z" in grid.xgcm_grid.axes: + axis_coords = grid.xgcm_grid.axes["Z"].coords.keys() + offsets["Z"] = 1 if "right" in axis_coords else 0 + else: + offsets["Z"] = 0 + return offsets + + def XLinear( particle_positions: dict[str, float | np.ndarray], grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]], @@ -165,14 +178,10 @@ def CGrid_Velocity( U = vectorfield.U.data V = vectorfield.V.data grid = vectorfield.grid + offsets = _get_offsets_dictionary(grid) tdim, zdim, ydim, xdim = U.shape[0], U.shape[1], U.shape[2], U.shape[3] lenT = 2 if np.any(tau > 0) else 1 - offsets = {} - for axis in ["X", "Y"]: - axis_coords = grid.xgcm_grid.axes[axis].coords.keys() - offsets[axis] = 1 if "right" in axis_coords else 0 - if grid.lon.ndim == 1: px = np.array([grid.lon[xi], grid.lon[xi + 1], grid.lon[xi + 1], grid.lon[xi]]) py = np.array([grid.lat[yi], grid.lat[yi], grid.lat[yi + 1], grid.lat[yi + 1]]) @@ -217,8 +226,8 @@ def _create_selection_dict(dims, zdir=False): if "Z" in axis_dim: if zdir: # Z coordinates: 1 point at zi and 1 point at zi+1 repeated for lenT time levels - zi_0 = np.clip(zi, 0, zdim - 1) - zi_1 = np.clip(zi + 1, 0, zdim - 1) + zi_0 = np.clip(zi + offsets["Z"], 0, zdim - 1) + zi_1 = np.clip(zi + offsets["Z"] + 1, 0, zdim - 1) zi_full = np.tile(np.array([zi_0, zi_1]).flatten(), lenT) else: # Z coordinates: 2 points at zi, repeated for lenT time levels @@ -335,6 +344,11 @@ def CGrid_Tracer( axis_dim = field.grid.get_axis_dim_mapping(field.data.dims) data = field.data + offsets = _get_offsets_dictionary(field.grid) + zi = np.clip(zi + offsets["Z"], 0, data.shape[1] - 1) + yi = np.clip(yi + offsets["Y"], 0, data.shape[2] - 1) + xi = np.clip(xi + offsets["X"], 0, data.shape[3] - 1) + lenT = 2 if np.any(tau > 0) else 1 if lenT == 2: From 24de000d4987cefd6899eb0df78a48dad678dfd1 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 7 Jan 2026 16:08:15 +0100 Subject: [PATCH 16/38] Move from_nemo body into parcels.convert --- src/parcels/_core/fieldset.py | 182 +----------------------------- src/parcels/convert.py | 207 ++++++++++++++++++++++++++++++++++ 2 files changed, 210 insertions(+), 179 deletions(-) create mode 100644 src/parcels/convert.py diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 2892779a1..90973235e 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -10,6 +10,7 @@ import xarray as xr import xgcm +from parcels import convert from parcels._core.converters import GeographicPolar from parcels._core.field import Field, VectorField from parcels._core.utils import sgrid @@ -20,6 +21,7 @@ from parcels._core.xgrid import _DEFAULT_XGCM_KWARGS, XGrid from parcels._logger import logger from parcels._typing import Mesh +from parcels.convert import _discover_U_and_V, _ds_rename_using_standard_names, _maybe_rename_coords from parcels.interpolators import ( CGrid_Velocity, UxPiecewiseConstantFace, @@ -267,47 +269,7 @@ def from_nemo(ds: xr.Dataset): If you encounter issues with your specific NEMO dataset, please open an issue on the Parcels GitHub repository with details about your dataset. """ - ds = ds.copy() - ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) - ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) - ds = _maybe_create_depth_dim(ds) - ds = _maybe_bring_UV_depths_to_depth(ds) - ds = _drop_unused_dimensions_and_coords(ds, _NEMO_DIMENSION_NAMES) - ds = _maybe_rename_coords(ds, _NEMO_AXIS_VARNAMES) - ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_NAMES) - ds = _set_coords(ds, _NEMO_DIMENSION_NAMES) - ds = _maybe_remove_depth_from_lonlat(ds) - ds = _set_axis_attrs(ds, _NEMO_AXIS_VARNAMES) - - expected_axes = set("XYZT") # TODO: Update after we have support for 2D spatial fields - if missing_axes := (expected_axes - set(ds.cf.axes)): - raise ValueError( - f"Dataset missing CF compliant metadata for axes " - f"{missing_axes}. Expected 'axis' attribute to be set " - f"on all dimension axes {expected_axes}. " - "HINT: Add xarray metadata attribute 'axis' to dimension - e.g., ds['lat'].attrs['axis'] = 'Y'" - ) - - if "W" in ds.data_vars: - # Negate W to convert from up positive to down positive (as that's the direction of positive z) - ds["W"].data *= -1 - if "grid" in ds.cf.cf_roles: - raise ValueError( - "Dataset already has a 'grid' variable (according to cf_roles). Didn't expect there to be grid metadata on copernicusmarine datasets - please open an issue with more information about your dataset." - ) - ds["grid"] = xr.DataArray( - 0, - attrs=sgrid.Grid2DMetadata( - cf_role="grid_topology", - topology_dimension=2, - node_dimensions=("glamf", "gphif"), - face_dimensions=( - sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), - sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), - ), - vertical_dimensions=(sgrid.DimDimPadding("z_center", "depth", sgrid.Padding.HIGH),), - ).to_attrs(), - ) + ds = convert.nemo_to_sgrid(ds) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") fieldset.V.units = GeographicPolar() if "UV" in fieldset.fields: @@ -513,134 +475,6 @@ def _format_calendar_error_message(field: Field, reference_datetime: TimeLike) - "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], } -_NEMO_CF_STANDARD_NAME_FALLBACKS = { - "UV": [ - ( - "sea_water_x_velocity", - "sea_water_y_velocity", - ), - ], - "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], -} - -_NEMO_DIMENSION_NAMES = ["x", "y", "time", "glamf", "gphif", "depth"] - -_NEMO_AXIS_VARNAMES = { - "X": "glamf", - "Y": "gphif", - "Z": "depth", - "T": "time", -} - -_NEMO_VARNAMES_MAPPING = { - "time_counter": "time", - "depthw": "depth", - "uo": "U", - "vo": "V", - "wo": "W", -} - - -def _maybe_bring_UV_depths_to_depth(ds): - if "U" in ds.variables and "depthu" in ds.U.coords and "depth" in ds.coords: - ds["U"] = ds["U"].assign_coords(depthu=ds["depth"].values).rename({"depthu": "depth"}) - if "V" in ds.variables and "depthv" in ds.V.coords and "depth" in ds.coords: - ds["V"] = ds["V"].assign_coords(depthv=ds["depth"].values).rename({"depthv": "depth"}) - return ds - - -def _maybe_create_depth_dim(ds): - if "depth" not in ds.dims: - ds = ds.expand_dims({"depth": [0]}) - ds["depth"] = xr.DataArray([0], dims=["depth"]) - return ds - - -def _maybe_rename_coords(ds, AXIS_VARNAMES): - try: - for axis, [coord] in ds.cf.axes.items(): - ds = ds.rename({coord: AXIS_VARNAMES[axis]}) - except ValueError as e: - raise ValueError(f"Multiple coordinates found on axis '{axis}'. Check your DataSet.") from e - return ds - - -def _maybe_rename_variables(ds, VARNAMES_MAPPING): - rename_dict = {old: new for old, new in VARNAMES_MAPPING.items() if (old in ds.data_vars) or (old in ds.coords)} - if rename_dict: - ds = ds.rename(rename_dict) - return ds - - -def _assign_dims_as_coords(ds, DIMENSION_NAMES): - for axis in DIMENSION_NAMES: - if axis in ds.dims and axis not in ds.coords: - ds = ds.assign_coords({axis: np.arange(ds.sizes[axis])}) - return ds - - -def _drop_unused_dimensions_and_coords(ds, DIMENSION_NAMES): - for dim in ds.dims: - if dim not in DIMENSION_NAMES: - ds = ds.drop_dims(dim, errors="ignore") - for coord in ds.coords: - if coord not in DIMENSION_NAMES: - ds = ds.drop_vars(coord, errors="ignore") - return ds - - -def _set_coords(ds, DIMENSION_NAMES): - for varname in DIMENSION_NAMES: - if varname in ds and varname not in ds.coords: - ds = ds.set_coords([varname]) - return ds - - -def _maybe_remove_depth_from_lonlat(ds): - for coord in ["glamf", "gphif"]: - if coord in ds.coords and "depth" in ds[coord].dims: - ds[coord] = ds[coord].squeeze("depth", drop=True) - return ds - - -def _set_axis_attrs(ds, AXIS_VARNAMES): - for axis, varname in AXIS_VARNAMES.items(): - if varname in ds.coords: - ds[varname].attrs["axis"] = axis - return ds - - -def _discover_U_and_V(ds: xr.Dataset, cf_standard_names_fallbacks) -> xr.Dataset: - # Assumes that the dataset has U and V data - - if "W" not in ds: - for cf_standard_name_W in cf_standard_names_fallbacks["W"]: - if cf_standard_name_W in ds.cf.standard_names: - ds = _ds_rename_using_standard_names(ds, {cf_standard_name_W: "W"}) - break - - if "U" in ds and "V" in ds: - return ds # U and V already present - elif "U" in ds or "V" in ds: - raise ValueError( - "Dataset has only one of the two variables 'U' and 'V'. Please rename the appropriate variable in your dataset to have both 'U' and 'V' for Parcels simulation." - ) - - for cf_standard_name_U, cf_standard_name_V in cf_standard_names_fallbacks["UV"]: - if cf_standard_name_U in ds.cf.standard_names: - if cf_standard_name_V not in ds.cf.standard_names: - raise ValueError( - f"Dataset has variable with CF standard name {cf_standard_name_U!r}, " - f"but not the matching variable with CF standard name {cf_standard_name_V!r}. " - "Please rename the appropriate variables in your dataset to have both 'U' and 'V' for Parcels simulation." - ) - else: - continue - - ds = _ds_rename_using_standard_names(ds, {cf_standard_name_U: "U", cf_standard_name_V: "V"}) - break - return ds - def _discover_fesom2_U_and_V(ds: ux.UxDataset) -> ux.UxDataset: # Common variable names for U and V found in UxDatasets @@ -684,16 +518,6 @@ def _discover_fesom2_U_and_V(ds: ux.UxDataset) -> ux.UxDataset: return ds -def _ds_rename_using_standard_names(ds: xr.Dataset | ux.UxDataset, name_dict: dict[str, str]) -> xr.Dataset: - for standard_name, rename_to in name_dict.items(): - name = ds.cf[standard_name].name - ds = ds.rename({name: rename_to}) - logger.info( - f"cf_xarray found variable {name!r} with CF standard name {standard_name!r} in dataset, renamed it to {rename_to!r} for Parcels simulation." - ) - return ds - - def _select_uxinterpolator(da: ux.UxDataArray): """Selects the appropriate uxarray interpolator for a given uxarray UxDataArray""" supported_uxinterp_mapping = { diff --git a/src/parcels/convert.py b/src/parcels/convert.py new file mode 100644 index 000000000..b57ad5263 --- /dev/null +++ b/src/parcels/convert.py @@ -0,0 +1,207 @@ +""" +This module provide a series of functions model outputs (which might be following their +own conventions) to metadata rich data (i.e., data following SGRID/UGRID as well as CF +conventions). Parcels needs this metadata rich data to discover grid geometries among other +things. + +These functions use knowledge about the model to attach any missing metadata. The functions +emit verbose messaging so that the user is kept in the loop. The returned output is an +Xarray dataset so that users can further provide any missing metadata that was unable to +be determined. +""" + +from __future__ import annotations + +import typing + +import numpy as np +import xarray as xr + +from parcels._core.utils import sgrid +from parcels._logger import logger + +if typing.TYPE_CHECKING: + import uxarray as ux + +_NEMO_CF_STANDARD_NAME_FALLBACKS = { + "UV": [ + ( + "sea_water_x_velocity", + "sea_water_y_velocity", + ), + ], + "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], +} + +_NEMO_DIMENSION_NAMES = ["x", "y", "time", "glamf", "gphif", "depth"] + +_NEMO_AXIS_VARNAMES = { + "X": "glamf", + "Y": "gphif", + "Z": "depth", + "T": "time", +} + +_NEMO_VARNAMES_MAPPING = { + "time_counter": "time", + "depthw": "depth", + "uo": "U", + "vo": "V", + "wo": "W", +} + + +def _maybe_bring_UV_depths_to_depth(ds): + if "U" in ds.variables and "depthu" in ds.U.coords and "depth" in ds.coords: + ds["U"] = ds["U"].assign_coords(depthu=ds["depth"].values).rename({"depthu": "depth"}) + if "V" in ds.variables and "depthv" in ds.V.coords and "depth" in ds.coords: + ds["V"] = ds["V"].assign_coords(depthv=ds["depth"].values).rename({"depthv": "depth"}) + return ds + + +def _maybe_create_depth_dim(ds): + if "depth" not in ds.dims: + ds = ds.expand_dims({"depth": [0]}) + ds["depth"] = xr.DataArray([0], dims=["depth"]) + return ds + + +def _maybe_rename_coords(ds, AXIS_VARNAMES): + try: + for axis, [coord] in ds.cf.axes.items(): + ds = ds.rename({coord: AXIS_VARNAMES[axis]}) + except ValueError as e: + raise ValueError(f"Multiple coordinates found on axis '{axis}'. Check your DataSet.") from e + return ds + + +def _maybe_rename_variables(ds, VARNAMES_MAPPING): + rename_dict = {old: new for old, new in VARNAMES_MAPPING.items() if (old in ds.data_vars) or (old in ds.coords)} + if rename_dict: + ds = ds.rename(rename_dict) + return ds + + +def _assign_dims_as_coords(ds, DIMENSION_NAMES): + for axis in DIMENSION_NAMES: + if axis in ds.dims and axis not in ds.coords: + ds = ds.assign_coords({axis: np.arange(ds.sizes[axis])}) + return ds + + +def _drop_unused_dimensions_and_coords(ds, DIMENSION_NAMES): + for dim in ds.dims: + if dim not in DIMENSION_NAMES: + ds = ds.drop_dims(dim, errors="ignore") + for coord in ds.coords: + if coord not in DIMENSION_NAMES: + ds = ds.drop_vars(coord, errors="ignore") + return ds + + +def _set_coords(ds, DIMENSION_NAMES): + for varname in DIMENSION_NAMES: + if varname in ds and varname not in ds.coords: + ds = ds.set_coords([varname]) + return ds + + +def _maybe_remove_depth_from_lonlat(ds): + for coord in ["glamf", "gphif"]: + if coord in ds.coords and "depth" in ds[coord].dims: + ds[coord] = ds[coord].squeeze("depth", drop=True) + return ds + + +def _set_axis_attrs(ds, AXIS_VARNAMES): + for axis, varname in AXIS_VARNAMES.items(): + if varname in ds.coords: + ds[varname].attrs["axis"] = axis + return ds + + +def _ds_rename_using_standard_names(ds: xr.Dataset | ux.UxDataset, name_dict: dict[str, str]) -> xr.Dataset: + for standard_name, rename_to in name_dict.items(): + name = ds.cf[standard_name].name + ds = ds.rename({name: rename_to}) + logger.info( + f"cf_xarray found variable {name!r} with CF standard name {standard_name!r} in dataset, renamed it to {rename_to!r} for Parcels simulation." + ) + return ds + + +def _discover_U_and_V(ds: xr.Dataset, cf_standard_names_fallbacks) -> xr.Dataset: + # Assumes that the dataset has U and V data + + if "W" not in ds: + for cf_standard_name_W in cf_standard_names_fallbacks["W"]: + if cf_standard_name_W in ds.cf.standard_names: + ds = _ds_rename_using_standard_names(ds, {cf_standard_name_W: "W"}) + break + + if "U" in ds and "V" in ds: + return ds # U and V already present + elif "U" in ds or "V" in ds: + raise ValueError( + "Dataset has only one of the two variables 'U' and 'V'. Please rename the appropriate variable in your dataset to have both 'U' and 'V' for Parcels simulation." + ) + + for cf_standard_name_U, cf_standard_name_V in cf_standard_names_fallbacks["UV"]: + if cf_standard_name_U in ds.cf.standard_names: + if cf_standard_name_V not in ds.cf.standard_names: + raise ValueError( + f"Dataset has variable with CF standard name {cf_standard_name_U!r}, " + f"but not the matching variable with CF standard name {cf_standard_name_V!r}. " + "Please rename the appropriate variables in your dataset to have both 'U' and 'V' for Parcels simulation." + ) + else: + continue + + ds = _ds_rename_using_standard_names(ds, {cf_standard_name_U: "U", cf_standard_name_V: "V"}) + break + return ds + + +def nemo_to_sgrid(ds: xr.Dataset): + ds = ds.copy() + ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) + ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) + ds = _maybe_create_depth_dim(ds) + ds = _maybe_bring_UV_depths_to_depth(ds) + ds = _drop_unused_dimensions_and_coords(ds, _NEMO_DIMENSION_NAMES) + ds = _maybe_rename_coords(ds, _NEMO_AXIS_VARNAMES) + ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_NAMES) + ds = _set_coords(ds, _NEMO_DIMENSION_NAMES) + ds = _maybe_remove_depth_from_lonlat(ds) + ds = _set_axis_attrs(ds, _NEMO_AXIS_VARNAMES) + + expected_axes = set("XYZT") # TODO: Update after we have support for 2D spatial fields + if missing_axes := (expected_axes - set(ds.cf.axes)): + raise ValueError( + f"Dataset missing CF compliant metadata for axes " + f"{missing_axes}. Expected 'axis' attribute to be set " + f"on all dimension axes {expected_axes}. " + "HINT: Add xarray metadata attribute 'axis' to dimension - e.g., ds['lat'].attrs['axis'] = 'Y'" + ) + + if "W" in ds.data_vars: + # Negate W to convert from up positive to down positive (as that's the direction of positive z) + ds["W"].data *= -1 + if "grid" in ds.cf.cf_roles: + raise ValueError( + "Dataset already has a 'grid' variable (according to cf_roles). Didn't expect there to be grid metadata on copernicusmarine datasets - please open an issue with more information about your dataset." + ) + ds["grid"] = xr.DataArray( + 0, + attrs=sgrid.Grid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("glamf", "gphif"), + face_dimensions=( + sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), + sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), + ), + vertical_dimensions=(sgrid.DimDimPadding("z_center", "depth", sgrid.Padding.HIGH),), + ).to_attrs(), + ) + return ds From 26c6d55d9d2a3a3a9c4ef732c5958b920d2bcfdd Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 7 Jan 2026 16:33:21 +0100 Subject: [PATCH 17/38] Move selection of vector interpolator to from_sgrid_conventions --- src/parcels/_core/fieldset.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 90973235e..f0e405f4d 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -272,10 +272,6 @@ def from_nemo(ds: xr.Dataset): ds = convert.nemo_to_sgrid(ds) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") fieldset.V.units = GeographicPolar() - if "UV" in fieldset.fields: - fieldset.UV.vector_interp_method = CGrid_Velocity - if "UVW" in fieldset.fields: - fieldset.UVW.vector_interp_method = CGrid_Velocity return fieldset def from_fesom2(ds: ux.UxDataset): @@ -395,11 +391,20 @@ def from_sgrid_conventions( if "U" in ds.data_vars and "V" in ds.data_vars: fields["U"] = Field("U", ds["U"], grid, XLinear) fields["V"] = Field("V", ds["V"], grid, XLinear) - fields["UV"] = VectorField("UV", fields["U"], fields["V"]) + fields["UV"] = VectorField( + "UV", + fields["U"], + fields["V"], + vector_interp_method=CGrid_Velocity, + # ^Seems to work with AGrid as well? (at least, tests aren't failing - + # either logic needs to be added to choose interpolator, or this interpolator should be renamed) + ) if "W" in ds.data_vars: fields["W"] = Field("W", ds["W"], grid, XLinear) - fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"]) + fields["UVW"] = VectorField( + "UVW", fields["U"], fields["V"], fields["W"], vector_interp_method=CGrid_Velocity + ) for varname in set(ds.data_vars) - set(fields.keys()) - skip_vars: fields[varname] = Field(varname, ds[varname], grid, XLinear) From c8800a2af0779c741b34a2f54f54b37b0ae29ff8 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 8 Jan 2026 13:06:53 +0100 Subject: [PATCH 18/38] Fix from_nemo --- src/parcels/_core/fieldset.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index f0e405f4d..b54cd1bd6 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -11,7 +11,7 @@ import xgcm from parcels import convert -from parcels._core.converters import GeographicPolar +from parcels._core.converters import Geographic, GeographicPolar from parcels._core.field import Field, VectorField from parcels._core.utils import sgrid from parcels._core.utils.string import _assert_str_and_python_varname @@ -409,7 +409,12 @@ def from_sgrid_conventions( for varname in set(ds.data_vars) - set(fields.keys()) - skip_vars: fields[varname] = Field(varname, ds[varname], grid, XLinear) - return FieldSet(list(fields.values())) + fieldset = FieldSet(list(fields.values())) + if mesh == "spherical": + fieldset.U.units = GeographicPolar() + fieldset.V.units = Geographic() + + return fieldset class CalendarError(Exception): # TODO: Move to a parcels errors module From ba7a95381957ed7fc7a92c36fa4f470aec6d0d77 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 8 Jan 2026 15:14:02 +0100 Subject: [PATCH 19/38] Refactor test_nemo_3D_curvilinear_fieldset --- tests/test_advection.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index 2a90cec4f..2ff9aa92d 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -24,7 +24,7 @@ AdvectionRK4_3D, AdvectionRK45, ) -from tests.utils import DEFAULT_PARTICLES, round_and_hash_float_array +from tests.utils import DEFAULT_PARTICLES @pytest.mark.parametrize("mesh", ["spherical", "flat"]) @@ -450,14 +450,18 @@ def test_nemo_3D_curvilinear_fieldset(kernel): fieldset = parcels.FieldSet.from_nemo(ds) npart = 10 - lons = np.linspace(1.9, 3.4, npart) - lats = np.linspace(52.5, 51.6, npart) - pset = parcels.ParticleSet(fieldset, lon=lons, lat=lats, z=np.ones_like(lons)) + lons_initial = np.linspace(1.9, 3.4, npart) + lats_initial = np.linspace(52.5, 51.6, npart) + z_initial = np.ones_like(lons_initial) + pset = parcels.ParticleSet(fieldset, lon=lons_initial, lat=lats_initial, z=z_initial) pset.execute(kernel, runtime=np.timedelta64(3, "D") + np.timedelta64(18, "h"), dt=np.timedelta64(6, "h")) if kernel == AdvectionRK4: - np.testing.assert_equal(round_and_hash_float_array([p.lon for p in pset], decimals=5), 29977383852960156017546) + np.testing.assert_allclose([p.z for p in pset], z_initial) elif kernel == AdvectionRK4_3D: # TODO check why decimals needs to be so low in RK4_3D (compare to v3) - np.testing.assert_equal(round_and_hash_float_array([p.z for p in pset], decimals=1), 29747210774230389239432) + np.testing.assert_allclose( + [p.z for p in pset], + [0.666162, 0.8667131, 0.92150104, 0.9605109, 0.9577529, 1.0041442, 1.0284728, 1.0033542, 1.2949713, 1.3928112], + ) # fmt:skip From 914e5ee22b42261d7a1ae73da946d3d94fb7fcd3 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 8 Jan 2026 15:25:51 +0100 Subject: [PATCH 20/38] Remove setting of fieldset.V.units in from_nemo It's not clear why this is here, nor why removing it causes a test failure. To be investigated another time. EDIT: This was introduced in https://github.com/Parcels-code/Parcels/commit/c311fba9a7ac3c204b54a5a51e3a04646263b692 - though we're investigating if this can be implemented another way since there should be no particular difference with NEMO --- src/parcels/_core/fieldset.py | 1 - tests/test_advection.py | 11 ++++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index b54cd1bd6..52a5800a6 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -271,7 +271,6 @@ def from_nemo(ds: xr.Dataset): """ ds = convert.nemo_to_sgrid(ds) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") - fieldset.V.units = GeographicPolar() return fieldset def from_fesom2(ds: ux.UxDataset): diff --git a/tests/test_advection.py b/tests/test_advection.py index 2ff9aa92d..ad325f74a 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -435,7 +435,16 @@ def test_nemo_curvilinear_fieldset(): np.testing.assert_allclose(pset.lat, latp, atol=1e-1) -@pytest.mark.parametrize("kernel", [AdvectionRK4, AdvectionRK4_3D]) +@pytest.mark.parametrize( + "kernel", + [ + AdvectionRK4, + pytest.mark.xfail( + AdvectionRK4_3D, + reason="from_nemo had 'fieldset.V.units = GeographicPolar()', I'm not sure _why_ this code is needed to get this to pass. To be further investigated.", + ), + ], +) def test_nemo_3D_curvilinear_fieldset(kernel): data_folder = parcels.download_example_dataset("NemoNorthSeaORCA025-N006_data") ds_fields = xr.open_mfdataset( From 5dc0a5f7f65ec21b9b2c4c7092da992e89645920 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 8 Jan 2026 16:25:05 +0100 Subject: [PATCH 21/38] Add discovery of `mesh` from convention metadata Gradually reducing the dependency on the `mesh` param --- src/parcels/_core/fieldset.py | 40 +++++++++++++++++++++++++++++++++-- src/parcels/convert.py | 4 ++++ 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 52a5800a6..801edb49b 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -270,7 +270,7 @@ def from_nemo(ds: xr.Dataset): """ ds = convert.nemo_to_sgrid(ds) - fieldset = FieldSet.from_sgrid_conventions(ds, mesh="spherical") + fieldset = FieldSet.from_sgrid_conventions(ds) return fieldset def from_fesom2(ds: ux.UxDataset): @@ -312,7 +312,7 @@ def from_fesom2(ds: ux.UxDataset): return FieldSet(list(fields.values())) def from_sgrid_conventions( - ds: xr.Dataset, mesh: Mesh + ds: xr.Dataset, mesh: Mesh | None = None ): # TODO: Update mesh to be discovered from the dataset metadata """Create a FieldSet from a dataset using SGRID convention metadata. @@ -342,6 +342,8 @@ def from_sgrid_conventions( See https://sgrid.github.io/ for more information on the SGRID conventions. """ ds = ds.copy() + if mesh is None: + mesh = _get_mesh_type_from_sgrid_dataset(ds) # Ensure time dimension has axis attribute if present if "time" in ds.dims and "time" in ds.coords: @@ -558,3 +560,37 @@ def _select_uxinterpolator(da: ux.UxDataArray): return supported_uxinterp_mapping[key] return None + + +# TODO: Refactor later into something like `parcels._metadata.discover(dataset)` helper that can be used to discover important metadata like this. I think this whole metadata handling should be refactored into its own module. +def _get_mesh_type_from_sgrid_dataset(ds_sgrid: xr.Dataset) -> Mesh: + """Small helper to inspect SGRID metadata and dataset metadata to determine mesh type.""" + grid_da = sgrid.get_grid_topology(ds_sgrid) + if grid_da is None: + raise ValueError("Dataset does not contain SGRID grid topology metadata (cf_role='grid_topology').") + + sgrid_metadata = sgrid.parse_grid_attrs(grid_da.attrs) + + fpoint_x, fpoint_y = sgrid_metadata.node_dimensions + + if _is_coordinate_in_degrees(ds_sgrid[fpoint_x]) ^ _is_coordinate_in_degrees(ds_sgrid[fpoint_x]): + msg = ( + f"Mismatch in units between X and Y coordinates.\n" + f" Coordinate {ds_sgrid[fpoint_x]!r} attrs: {ds_sgrid[fpoint_x].attrs}\n" + f" Coordinate {ds_sgrid[fpoint_y]!r} attrs: {ds_sgrid[fpoint_y].attrs}\n" + ) + raise ValueError(msg) + + return "spherical" if _is_coordinate_in_degrees(ds_sgrid[fpoint_x]) else "flat" + + +def _is_coordinate_in_degrees(da: xr.DataArray) -> bool: + match da.attrs.get("units"): + case None: + raise ValueError( + f"Coordinate {da.name} of your dataset has no 'units' attribute - we don't know what the spatial units are." + ) + case "degrees": + return True + case _: + return False diff --git a/src/parcels/convert.py b/src/parcels/convert.py index b57ad5263..d22c7a9f4 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -204,4 +204,8 @@ def nemo_to_sgrid(ds: xr.Dataset): vertical_dimensions=(sgrid.DimDimPadding("z_center", "depth", sgrid.Padding.HIGH),), ).to_attrs(), ) + + # NEMO models are always in degrees + ds["glamf"].attrs["units"] = "degrees" + ds["gphif"].attrs["units"] = "degrees" return ds From ea4e99358e511bbc4035149e4ac6019e9c664ada Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 8 Jan 2026 16:33:34 +0100 Subject: [PATCH 22/38] Fix xfail mark --- tests/test_advection.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index ad325f74a..0c490958c 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -439,9 +439,11 @@ def test_nemo_curvilinear_fieldset(): "kernel", [ AdvectionRK4, - pytest.mark.xfail( + pytest.param( AdvectionRK4_3D, - reason="from_nemo had 'fieldset.V.units = GeographicPolar()', I'm not sure _why_ this code is needed to get this to pass. To be further investigated.", + marks=pytest.mark.xfail( + reason="from_nemo had 'fieldset.V.units = GeographicPolar()', I'm not sure _why_ this code is needed to get this to pass. To be further investigated.", + ), ), ], ) From fd699f279ba00ad6c30c02cee4ef48e4e900e9e1 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 8 Jan 2026 17:29:34 +0100 Subject: [PATCH 23/38] Add helpers get_value_by_id --- src/parcels/_core/utils/sgrid.py | 50 ++++++++++++++++++++++++++++++++ tests/utils/test_sgrid.py | 24 +++++++++++++++ 2 files changed, 74 insertions(+) diff --git a/src/parcels/_core/utils/sgrid.py b/src/parcels/_core/utils/sgrid.py index ff87d5081..5ecd7cec4 100644 --- a/src/parcels/_core/utils/sgrid.py +++ b/src/parcels/_core/utils/sgrid.py @@ -47,6 +47,30 @@ def to_attrs(self) -> dict[str, str | int]: ... def from_attrs(cls, d: dict[str, Hashable]) -> Self: ... +# Note that - for some optional attributes in the SGRID spec - these IDs are not available +# hence this isn't full coverage +_ID_FETCHERS_GRID2DMETADATA = { + "node_dimension1": lambda meta: meta.node_dimensions[0], + "node_dimension2": lambda meta: meta.node_dimensions[1], + "face_dimension1": lambda meta: meta.face_dimensions[0].dim1, + "face_dimension2": lambda meta: meta.face_dimensions[1].dim1, + "type1": lambda meta: meta.face_dimensions[0].padding, + "type2": lambda meta: meta.face_dimensions[1].padding, +} + +_ID_FETCHERS_GRID3DMETADATA = { + "node_dimension1": lambda meta: meta.node_dimensions[0], + "node_dimension2": lambda meta: meta.node_dimensions[1], + "node_dimension3": lambda meta: meta.node_dimensions[2], + "face_dimension1": lambda meta: meta.volume_dimensions[0].dim1, + "face_dimension2": lambda meta: meta.volume_dimensions[1].dim1, + "face_dimension3": lambda meta: meta.volume_dimensions[2].dim1, + "type1": lambda meta: meta.volume_dimensions[0].padding, + "type2": lambda meta: meta.volume_dimensions[1].padding, + "type3": lambda meta: meta.volume_dimensions[2].padding, +} + + class Grid2DMetadata(AttrsSerializable): def __init__( self, @@ -152,6 +176,19 @@ def to_attrs(self) -> dict[str, str | int]: def rename_dims(self, dims_dict: dict[str, str]) -> Self: return _metadata_rename_dims(self, dims_dict) + def get_value_by_id(self, id: str) -> str: + """In the SGRID specification different parts of the spec are identified to with different "ID"s. + + Easily extract the value for a given ID. + + Example + ------- + # Get padding 2 + >>> get_name_from_id("type2") + "low" + """ + return _ID_FETCHERS_GRID2DMETADATA[id](self) + class Grid3DMetadata(AttrsSerializable): def __init__( @@ -249,6 +286,19 @@ def to_attrs(self) -> dict[str, str | int]: def rename_dims(self, dims_dict: dict[str, str]) -> Self: return _metadata_rename_dims(self, dims_dict) + def get_value_by_id(self, id: str) -> str: + """In the SGRID specification different parts of the spec are identified to with different "ID"s. + + Easily extract the value for a given ID. + + Example + ------- + # Get padding 2 + >>> get_name_from_id("type2") + "low" + """ + return _ID_FETCHERS_GRID3DMETADATA[id](self) + @dataclass class DimDimPadding: diff --git a/tests/utils/test_sgrid.py b/tests/utils/test_sgrid.py index e73c97865..59c5f6fee 100644 --- a/tests/utils/test_sgrid.py +++ b/tests/utils/test_sgrid.py @@ -52,6 +52,30 @@ def create_example_grid3dmetadata(with_node_coordinates: bool): grid3dmetadata = create_example_grid3dmetadata(with_node_coordinates=True) +@pytest.mark.parametrize( + ("sgrid_metadata", "id", "value"), + [ + (grid2dmetadata, "node_dimension1", "node_dimension1"), + (grid2dmetadata, "node_dimension2", "node_dimension2"), + (grid2dmetadata, "face_dimension1", "face_dimension1"), + (grid2dmetadata, "face_dimension2", "face_dimension2"), + # (grid2dmetadata, "vertical_dimensions_dim1", "vertical_dimensions_dim1"), #! ID's NOT IMPLEMENTED IN SGRID SPEC + # (grid2dmetadata, "vertical_dimensions_dim2", "vertical_dimensions_dim2"), + (grid3dmetadata, "node_dimension1", "node_dimension1"), + (grid3dmetadata, "node_dimension2", "node_dimension2"), + (grid3dmetadata, "node_dimension3", "node_dimension3"), + (grid3dmetadata, "face_dimension1", "face_dimension1"), + (grid3dmetadata, "face_dimension2", "face_dimension2"), + (grid3dmetadata, "face_dimension3", "face_dimension3"), + (grid3dmetadata, "type1", sgrid.Padding.LOW), + (grid3dmetadata, "type2", sgrid.Padding.LOW), + (grid3dmetadata, "type3", sgrid.Padding.LOW), + ], +) +def test_get_value_by_id(sgrid_metadata: sgrid.Grid2DMetadata | sgrid.Grid3DMetadata, id, value): + assert sgrid_metadata.get_value_by_id(id) == value + + def dummy_sgrid_ds(grid: sgrid.Grid2DMetadata | sgrid.Grid3DMetadata) -> xr.Dataset: if isinstance(grid, sgrid.Grid2DMetadata): return dummy_sgrid_2d_ds(grid) From 875fe702e1c3b09d82479b6eae6795991bdaad5f Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 8 Jan 2026 17:51:51 +0100 Subject: [PATCH 24/38] Add test_convert.test_nemo_to_sgrid --- src/parcels/convert.py | 3 ++- tests/test_convert.py | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 1 deletion(-) create mode 100644 tests/test_convert.py diff --git a/src/parcels/convert.py b/src/parcels/convert.py index d22c7a9f4..4ab073c0f 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -196,7 +196,8 @@ def nemo_to_sgrid(ds: xr.Dataset): attrs=sgrid.Grid2DMetadata( cf_role="grid_topology", topology_dimension=2, - node_dimensions=("glamf", "gphif"), + node_dimensions=("x", "y"), + node_coordinates=("glamf", "gphif"), face_dimensions=( sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), diff --git a/tests/test_convert.py b/tests/test_convert.py new file mode 100644 index 000000000..e0e6008af --- /dev/null +++ b/tests/test_convert.py @@ -0,0 +1,39 @@ +import xarray as xr + +import parcels +import parcels.convert as convert +from parcels._core.utils import sgrid + + +def test_nemo_to_sgrid(): + data_folder = parcels.download_example_dataset("NemoCurvilinear_data") + ds_fields = xr.open_mfdataset( + data_folder.glob("*.nc4"), + data_vars="minimal", + coords="minimal", + compat="override", + ) + ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True) + ds = convert.nemo_to_sgrid(ds_fields) + + assert ds["grid"].attrs == { + "cf_role": "grid_topology", + "topology_dimension": 2, + "node_dimensions": "x y", + "face_dimensions": "x_center:x (padding:low) y_center:y (padding:low)", + "node_coordinates": "glamf gphif", + "vertical_dimensions": "z_center:depth (padding:high)", + } + + meta = sgrid.parse_grid_attrs(ds["grid"].attrs) + + # Assuming that node_dimension1 and node_dimension2 correspond to X and Y respectively + # check that U and V are properly defined on the staggered grid + assert { + meta.get_value_by_id("node_dimension1"), # X edge + meta.get_value_by_id("face_dimension2"), # Y center + } in set(ds["U"].dims) + assert { + meta.get_value_by_id("face_dimension1"), # X center + meta.get_value_by_id("node_dimension2"), # Y edge + } in set(ds["V"].dims) From 31c63dd3653f2322ecba9350a648c373940f3e66 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 9 Jan 2026 09:58:57 +0100 Subject: [PATCH 25/38] Fix param naming to obey pep8 --- src/parcels/convert.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/parcels/convert.py b/src/parcels/convert.py index 4ab073c0f..be657e746 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -66,41 +66,41 @@ def _maybe_create_depth_dim(ds): return ds -def _maybe_rename_coords(ds, AXIS_VARNAMES): +def _maybe_rename_coords(ds, axis_varnames): try: for axis, [coord] in ds.cf.axes.items(): - ds = ds.rename({coord: AXIS_VARNAMES[axis]}) + ds = ds.rename({coord: axis_varnames[axis]}) except ValueError as e: raise ValueError(f"Multiple coordinates found on axis '{axis}'. Check your DataSet.") from e return ds -def _maybe_rename_variables(ds, VARNAMES_MAPPING): - rename_dict = {old: new for old, new in VARNAMES_MAPPING.items() if (old in ds.data_vars) or (old in ds.coords)} +def _maybe_rename_variables(ds, varnames_mapping): + rename_dict = {old: new for old, new in varnames_mapping.items() if (old in ds.data_vars) or (old in ds.coords)} if rename_dict: ds = ds.rename(rename_dict) return ds -def _assign_dims_as_coords(ds, DIMENSION_NAMES): - for axis in DIMENSION_NAMES: +def _assign_dims_as_coords(ds, dimension_names): + for axis in dimension_names: if axis in ds.dims and axis not in ds.coords: ds = ds.assign_coords({axis: np.arange(ds.sizes[axis])}) return ds -def _drop_unused_dimensions_and_coords(ds, DIMENSION_NAMES): +def _drop_unused_dimensions_and_coords(ds, dimension_names): for dim in ds.dims: - if dim not in DIMENSION_NAMES: + if dim not in dimension_names: ds = ds.drop_dims(dim, errors="ignore") for coord in ds.coords: - if coord not in DIMENSION_NAMES: + if coord not in dimension_names: ds = ds.drop_vars(coord, errors="ignore") return ds -def _set_coords(ds, DIMENSION_NAMES): - for varname in DIMENSION_NAMES: +def _set_coords(ds, dimension_names): + for varname in dimension_names: if varname in ds and varname not in ds.coords: ds = ds.set_coords([varname]) return ds @@ -113,8 +113,8 @@ def _maybe_remove_depth_from_lonlat(ds): return ds -def _set_axis_attrs(ds, AXIS_VARNAMES): - for axis, varname in AXIS_VARNAMES.items(): +def _set_axis_attrs(ds, axis_varnames): + for axis, varname in axis_varnames.items(): if varname in ds.coords: ds[varname].attrs["axis"] = axis return ds From 24357bfb556a44442f1ebf154ad393a65621bde5 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 9 Jan 2026 15:37:13 +0100 Subject: [PATCH 26/38] Update API for convert.nemo_to_sgrid Updates the API for conversion to be more closely aligned with the input data. Also handles the U and V fields separately - correctly assigning the dimension naming before merging into a single dataset. --- src/parcels/convert.py | 66 ++++++++++++++++++++++++++++++------------ tests/test_convert.py | 17 +++++------ 2 files changed, 55 insertions(+), 28 deletions(-) diff --git a/src/parcels/convert.py b/src/parcels/convert.py index be657e746..e9313fb03 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -33,13 +33,15 @@ "W": ["upward_sea_water_velocity", "vertical_sea_water_velocity"], } -_NEMO_DIMENSION_NAMES = ["x", "y", "time", "glamf", "gphif", "depth"] +_NEMO_DIMENSION_COORD_NAMES = ["x", "y", "time", "x", "x_center", "y", "y_center", "depth", "glamf", "gphif"] _NEMO_AXIS_VARNAMES = { - "X": "glamf", - "Y": "gphif", - "Z": "depth", - "T": "time", + "x": "X", + "x_center": "X", + "y": "Y", + "y_center": "Y", + "depth": "Z", + "time": "T", } _NEMO_VARNAMES_MAPPING = { @@ -89,12 +91,12 @@ def _assign_dims_as_coords(ds, dimension_names): return ds -def _drop_unused_dimensions_and_coords(ds, dimension_names): +def _drop_unused_dimensions_and_coords(ds, dimension_and_coord_names): for dim in ds.dims: - if dim not in dimension_names: + if dim not in dimension_and_coord_names: ds = ds.drop_dims(dim, errors="ignore") for coord in ds.coords: - if coord not in dimension_names: + if coord not in dimension_and_coord_names: ds = ds.drop_vars(coord, errors="ignore") return ds @@ -113,10 +115,9 @@ def _maybe_remove_depth_from_lonlat(ds): return ds -def _set_axis_attrs(ds, axis_varnames): - for axis, varname in axis_varnames.items(): - if varname in ds.coords: - ds[varname].attrs["axis"] = axis +def _set_axis_attrs(ds, dim_axis): + for dim, axis in dim_axis.items(): + ds[dim].attrs["axis"] = axis return ds @@ -162,16 +163,45 @@ def _discover_U_and_V(ds: xr.Dataset, cf_standard_names_fallbacks) -> xr.Dataset return ds -def nemo_to_sgrid(ds: xr.Dataset): - ds = ds.copy() +def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset]): + fields = fields.copy() + coords = coords[["gphif", "glamf"]] + + for name, field_da in fields.items(): + if isinstance(field_da, xr.Dataset): + field_da = field_da[name] + # TODO: logging message, warn if multiple fields are in this dataset + + match name: + case "U": + field_da = field_da.rename({"y": "y_center"}) + case "V": + field_da = field_da.rename({"x": "x_center"}) + case _: + pass + + fields[name] = field_da + + if "time" in coords.dims: + if coords.dims["time"] != 1: + raise ValueError("Time dimension in coords must be length 1 (i.e., no time-varying grid).") + coords = coords.isel(time=0).drop("time") + if len(coords.dims) == 3: + for dim, len_ in coords.dims.items(): + if len_ == 1: + # TODO: log statement about selecting along z dim of 1 + coords = coords.isel({dim: 0}) + if len(coords.dims) != 2: + raise ValueError("Expected coordsinates to be 2 dimensional") + + ds = xr.merge(list(fields.values()) + [coords]) ds = _maybe_rename_variables(ds, _NEMO_VARNAMES_MAPPING) ds = _discover_U_and_V(ds, _NEMO_CF_STANDARD_NAME_FALLBACKS) ds = _maybe_create_depth_dim(ds) ds = _maybe_bring_UV_depths_to_depth(ds) - ds = _drop_unused_dimensions_and_coords(ds, _NEMO_DIMENSION_NAMES) - ds = _maybe_rename_coords(ds, _NEMO_AXIS_VARNAMES) - ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_NAMES) - ds = _set_coords(ds, _NEMO_DIMENSION_NAMES) + ds = _drop_unused_dimensions_and_coords(ds, _NEMO_DIMENSION_COORD_NAMES) + ds = _assign_dims_as_coords(ds, _NEMO_DIMENSION_COORD_NAMES) + ds = _set_coords(ds, _NEMO_DIMENSION_COORD_NAMES) ds = _maybe_remove_depth_from_lonlat(ds) ds = _set_axis_attrs(ds, _NEMO_AXIS_VARNAMES) diff --git a/tests/test_convert.py b/tests/test_convert.py index e0e6008af..d7777d49c 100644 --- a/tests/test_convert.py +++ b/tests/test_convert.py @@ -7,14 +7,11 @@ def test_nemo_to_sgrid(): data_folder = parcels.download_example_dataset("NemoCurvilinear_data") - ds_fields = xr.open_mfdataset( - data_folder.glob("*.nc4"), - data_vars="minimal", - coords="minimal", - compat="override", - ) - ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True) - ds = convert.nemo_to_sgrid(ds_fields) + U = xr.open_mfdataset(data_folder.glob("*U.nc4")) + V = xr.open_mfdataset(data_folder.glob("*V.nc4")) + coords = xr.open_dataset(data_folder / "mesh_mask.nc4") + + ds = convert.nemo_to_sgrid(U=U, V=V, coords=coords) assert ds["grid"].attrs == { "cf_role": "grid_topology", @@ -32,8 +29,8 @@ def test_nemo_to_sgrid(): assert { meta.get_value_by_id("node_dimension1"), # X edge meta.get_value_by_id("face_dimension2"), # Y center - } in set(ds["U"].dims) + }.issubset(set(ds["U"].dims)) assert { meta.get_value_by_id("face_dimension1"), # X center meta.get_value_by_id("node_dimension2"), # Y edge - } in set(ds["V"].dims) + }.issubset(set(ds["V"].dims)) From 3670dd524305e74be6bb9c30ef349bc59d1f1a5d Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 9 Jan 2026 16:05:40 +0100 Subject: [PATCH 27/38] Remove FieldSet.from_nemo Update tests --- src/parcels/_core/fieldset.py | 30 ++---------------------------- src/parcels/convert.py | 27 +++++++++++++++++++++++++-- tests/test_advection.py | 33 +++++++++++++++------------------ 3 files changed, 42 insertions(+), 48 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 801edb49b..f36b97426 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -10,7 +10,6 @@ import xarray as xr import xgcm -from parcels import convert from parcels._core.converters import Geographic, GeographicPolar from parcels._core.field import Field, VectorField from parcels._core.utils import sgrid @@ -248,31 +247,6 @@ def from_copernicusmarine(ds: xr.Dataset): ) return FieldSet.from_sgrid_conventions(ds, mesh="spherical") - def from_nemo(ds: xr.Dataset): - """Create a FieldSet from a xarray.Dataset from NEMO netcdf files. - - Parameters - ---------- - ds : xarray.Dataset - xarray.Dataset as obtained from a set of NEMO netcdf files. - - Returns - ------- - FieldSet - FieldSet object containing the fields from the dataset that can be used for a Parcels simulation. - - Notes - ----- - The NEMO model (https://www.nemo-ocean.eu/) is used by a variety of oceanographic institutions around the world. - Output from these models may differ subtly in terms of variable names and metadata conventions. - This function attempts to standardize these differences to create a Parcels FieldSet. - If you encounter issues with your specific NEMO dataset, please open an issue on the Parcels GitHub repository with details about your dataset. - - """ - ds = convert.nemo_to_sgrid(ds) - fieldset = FieldSet.from_sgrid_conventions(ds) - return fieldset - def from_fesom2(ds: ux.UxDataset): """Create a FieldSet from a FESOM2 uxarray.UxDataset. @@ -571,7 +545,7 @@ def _get_mesh_type_from_sgrid_dataset(ds_sgrid: xr.Dataset) -> Mesh: sgrid_metadata = sgrid.parse_grid_attrs(grid_da.attrs) - fpoint_x, fpoint_y = sgrid_metadata.node_dimensions + fpoint_x, fpoint_y = sgrid_metadata.node_coordinates if _is_coordinate_in_degrees(ds_sgrid[fpoint_x]) ^ _is_coordinate_in_degrees(ds_sgrid[fpoint_x]): msg = ( @@ -588,7 +562,7 @@ def _is_coordinate_in_degrees(da: xr.DataArray) -> bool: match da.attrs.get("units"): case None: raise ValueError( - f"Coordinate {da.name} of your dataset has no 'units' attribute - we don't know what the spatial units are." + f"Coordinate {da.name!r} of your dataset has no 'units' attribute - we don't know what the spatial units are." ) case "degrees": return True diff --git a/src/parcels/convert.py b/src/parcels/convert.py index e9313fb03..bc3219ee0 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -163,7 +163,28 @@ def _discover_U_and_V(ds: xr.Dataset, cf_standard_names_fallbacks) -> xr.Dataset return ds -def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset]): +def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset | xr.DataArray]): + # TODO: Update docstring + """Create a FieldSet from a xarray.Dataset from NEMO netcdf files. + + Parameters + ---------- + ds : xarray.Dataset + xarray.Dataset as obtained from a set of NEMO netcdf files. + + Returns + ------- + xarray.Dataset + Dataset object following SGRID conventions to be (optionally) modified and passed to a FieldSet constructor. + + Notes + ----- + The NEMO model (https://www.nemo-ocean.eu/) is used by a variety of oceanographic institutions around the world. + Output from these models may differ subtly in terms of variable names and metadata conventions. + This function attempts to standardize these differences to create a Parcels FieldSet. + If you encounter issues with your specific NEMO dataset, please open an issue on the Parcels GitHub repository with details about your dataset. + + """ fields = fields.copy() coords = coords[["gphif", "glamf"]] @@ -171,6 +192,8 @@ def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset]): if isinstance(field_da, xr.Dataset): field_da = field_da[name] # TODO: logging message, warn if multiple fields are in this dataset + else: + field_da = field_da.rename(name) match name: case "U": @@ -179,7 +202,7 @@ def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset]): field_da = field_da.rename({"x": "x_center"}) case _: pass - + field_da = field_da.reset_coords(drop=True) fields[name] = field_da if "time" in coords.dims: diff --git a/tests/test_advection.py b/tests/test_advection.py index 0c490958c..3592954f0 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -416,14 +416,13 @@ def UpdateP(particles, fieldset): # pragma: no cover def test_nemo_curvilinear_fieldset(): data_folder = parcels.download_example_dataset("NemoCurvilinear_data") - ds_fields = xr.open_mfdataset( - data_folder.glob("*.nc4"), - data_vars="minimal", - coords="minimal", - compat="override", - ) - ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True) - fieldset = parcels.FieldSet.from_nemo(ds_fields) + U = xr.open_mfdataset(data_folder.glob("*U.nc4")) + V = xr.open_mfdataset(data_folder.glob("*V.nc4")) + coords = xr.open_dataset(data_folder / "mesh_mask.nc4") + + ds = parcels.convert.nemo_to_sgrid(U=U, V=V, coords=coords) + + fieldset = parcels.FieldSet.from_sgrid_conventions(ds) npart = 20 lonp = 30 * np.ones(npart) @@ -449,16 +448,14 @@ def test_nemo_curvilinear_fieldset(): ) def test_nemo_3D_curvilinear_fieldset(kernel): data_folder = parcels.download_example_dataset("NemoNorthSeaORCA025-N006_data") - ds_fields = xr.open_mfdataset( - data_folder.glob("ORCA*.nc"), - data_vars="minimal", - coords="minimal", - compat="override", - ) - ds_coords = xr.open_dataset(data_folder / "coordinates.nc", decode_times=False) - ds_coords = ds_coords.isel(time=0, drop=True) - ds = xr.merge([ds_fields, ds_coords[["glamf", "gphif"]]]) - fieldset = parcels.FieldSet.from_nemo(ds) + U = xr.open_mfdataset(data_folder.glob("*U.nc")) + V = xr.open_mfdataset(data_folder.glob("*V.nc")) + W = xr.open_mfdataset(data_folder.glob("*W.nc")) + coords = xr.open_dataset(data_folder / "coordinates.nc", decode_times=False) + + ds = parcels.convert.nemo_to_sgrid(U=U["uo"], V=V["vo"], W=W["wo"], coords=coords) + + fieldset = parcels.FieldSet.from_sgrid_conventions(ds) npart = 10 lons_initial = np.linspace(1.9, 3.4, npart) From 0b9c787c8d14cf02cc393aff805142308c2c3398 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 9 Jan 2026 16:22:22 +0100 Subject: [PATCH 28/38] Fix from_sgrid_conventions Avoid assuming there's a U and V field. Maybe this should be refactored later... --- src/parcels/_core/fieldset.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index f36b97426..5b687090b 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -386,8 +386,10 @@ def from_sgrid_conventions( fieldset = FieldSet(list(fields.values())) if mesh == "spherical": - fieldset.U.units = GeographicPolar() - fieldset.V.units = Geographic() + if "U" in fieldset.fields: + fieldset.U.units = GeographicPolar() + if "V" in fieldset.fields: + fieldset.V.units = Geographic() return fieldset From 1926f6a9602d7ed65765a78180b6b6f538549006 Mon Sep 17 00:00:00 2001 From: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 12 Jan 2026 10:19:13 +0100 Subject: [PATCH 29/38] Update src/parcels/_core/utils/sgrid.py Co-authored-by: Erik van Sebille --- src/parcels/_core/utils/sgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/utils/sgrid.py b/src/parcels/_core/utils/sgrid.py index 5ecd7cec4..b229ba929 100644 --- a/src/parcels/_core/utils/sgrid.py +++ b/src/parcels/_core/utils/sgrid.py @@ -287,7 +287,7 @@ def rename_dims(self, dims_dict: dict[str, str]) -> Self: return _metadata_rename_dims(self, dims_dict) def get_value_by_id(self, id: str) -> str: - """In the SGRID specification different parts of the spec are identified to with different "ID"s. + """In the SGRID specification for 3D grids, different parts of the spec are identified by different "ID"s. Easily extract the value for a given ID. From 8f76e8555780d662237a666dc15a57f062cf906a Mon Sep 17 00:00:00 2001 From: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 12 Jan 2026 10:20:20 +0100 Subject: [PATCH 30/38] Update src/parcels/convert.py Co-authored-by: Erik van Sebille --- src/parcels/convert.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/convert.py b/src/parcels/convert.py index bc3219ee0..a04870981 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -7,7 +7,7 @@ These functions use knowledge about the model to attach any missing metadata. The functions emit verbose messaging so that the user is kept in the loop. The returned output is an Xarray dataset so that users can further provide any missing metadata that was unable to -be determined. +be determined before they pass it to the FieldSet constructor. """ from __future__ import annotations From 4ff5f109023130464d9b04f41bab4852a5267c17 Mon Sep 17 00:00:00 2001 From: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 12 Jan 2026 10:45:45 +0100 Subject: [PATCH 31/38] Update src/parcels/_core/utils/sgrid.py Co-authored-by: Erik van Sebille --- src/parcels/_core/utils/sgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/utils/sgrid.py b/src/parcels/_core/utils/sgrid.py index b229ba929..0fe3f26c8 100644 --- a/src/parcels/_core/utils/sgrid.py +++ b/src/parcels/_core/utils/sgrid.py @@ -177,7 +177,7 @@ def rename_dims(self, dims_dict: dict[str, str]) -> Self: return _metadata_rename_dims(self, dims_dict) def get_value_by_id(self, id: str) -> str: - """In the SGRID specification different parts of the spec are identified to with different "ID"s. + """In the SGRID specification for 2D grids, different parts of the spec are identified by different "ID"s. Easily extract the value for a given ID. From 6e133c91538ba6af7659d1bad204fb7f90aab2dc Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 12 Jan 2026 14:21:31 +0100 Subject: [PATCH 32/38] Update docs/user_guide/index.md --- docs/user_guide/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index cca43eee1..04d004884 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -30,7 +30,7 @@ TODO: Add links to Reference API throughout :titlesonly: examples/explanation_grids.md examples/tutorial_nemo_curvilinear.ipynb -examples/tutorial_nemo_3D.ipynb +examples_v3/tutorial_nemo_3D.ipynb # TODO move to examples folder examples/tutorial_unitconverters.ipynb examples/tutorial_nestedgrids.ipynb ``` From cd602977d1a7b426b418775fef42a05a15cb11c7 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 12 Jan 2026 14:39:03 +0100 Subject: [PATCH 33/38] Avoid xarray FutureWaening Avoiding a FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`. --- src/parcels/convert.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parcels/convert.py b/src/parcels/convert.py index a04870981..637f52700 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -206,11 +206,11 @@ def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset | xr.Dat fields[name] = field_da if "time" in coords.dims: - if coords.dims["time"] != 1: + if coords.sizes["time"] != 1: raise ValueError("Time dimension in coords must be length 1 (i.e., no time-varying grid).") coords = coords.isel(time=0).drop("time") if len(coords.dims) == 3: - for dim, len_ in coords.dims.items(): + for dim, len_ in coords.sizes.items(): if len_ == 1: # TODO: log statement about selecting along z dim of 1 coords = coords.isel({dim: 0}) From e912fd70d86187d01cf2d9350492675732e4086a Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 13 Jan 2026 13:00:49 +0100 Subject: [PATCH 34/38] Fix glamf and gphif --- .../examples/tutorial_nemo_curvilinear.ipynb | 12 +++++++----- src/parcels/convert.py | 10 +++++++--- tests/test_convert.py | 2 +- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index 3defa1f8a..24422013a 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb @@ -71,10 +71,12 @@ " compat=\"override\",\n", ")\n", "\n", - "# Drop some auxiliary dimensions - otherwise Parcels will complain\n", - "ds_fields = ds_fields.isel(time=0, z_a=0, z=0, drop=True)\n", + "ds_coords = xr.open_dataset(data_folder / \"mesh_mask.nc4\", decode_times=False)\n", + "ds_fset = parcels.convert.nemo_to_sgrid(\n", + " U=ds_fields[\"U\"], V=ds_fields[\"V\"], coords=ds_coords\n", + ")\n", "\n", - "fieldset = parcels.FieldSet.from_nemo(ds_fields)" + "fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)" ] }, { @@ -272,7 +274,7 @@ ], "metadata": { "kernelspec": { - "display_name": "default", + "display_name": "test-notebooks-latest", "language": "python", "name": "python3" }, @@ -286,7 +288,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.14.2" + "version": "3.13.9" } }, "nbformat": 4, diff --git a/src/parcels/convert.py b/src/parcels/convert.py index 637f52700..c38183621 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -244,13 +244,17 @@ def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset | xr.Dat raise ValueError( "Dataset already has a 'grid' variable (according to cf_roles). Didn't expect there to be grid metadata on copernicusmarine datasets - please open an issue with more information about your dataset." ) + + # Update to use lon and lat internally + + ds = ds.rename({"gphif": "lat", "glamf": "lon"}) # TODO: Logging message about rename ds["grid"] = xr.DataArray( 0, attrs=sgrid.Grid2DMetadata( cf_role="grid_topology", topology_dimension=2, node_dimensions=("x", "y"), - node_coordinates=("glamf", "gphif"), + node_coordinates=("lon", "lat"), face_dimensions=( sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), @@ -260,6 +264,6 @@ def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset | xr.Dat ) # NEMO models are always in degrees - ds["glamf"].attrs["units"] = "degrees" - ds["gphif"].attrs["units"] = "degrees" + ds["lon"].attrs["units"] = "degrees" + ds["lat"].attrs["units"] = "degrees" return ds diff --git a/tests/test_convert.py b/tests/test_convert.py index d7777d49c..11a7029dc 100644 --- a/tests/test_convert.py +++ b/tests/test_convert.py @@ -18,7 +18,7 @@ def test_nemo_to_sgrid(): "topology_dimension": 2, "node_dimensions": "x y", "face_dimensions": "x_center:x (padding:low) y_center:y (padding:low)", - "node_coordinates": "glamf gphif", + "node_coordinates": "lon lat", "vertical_dimensions": "z_center:depth (padding:high)", } From e380c0a324f55789062ae1e60647935a68447a2d Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 13 Jan 2026 13:06:57 +0100 Subject: [PATCH 35/38] Refactor --- src/parcels/convert.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/parcels/convert.py b/src/parcels/convert.py index c38183621..d2c7533dd 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -245,16 +245,13 @@ def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset | xr.Dat "Dataset already has a 'grid' variable (according to cf_roles). Didn't expect there to be grid metadata on copernicusmarine datasets - please open an issue with more information about your dataset." ) - # Update to use lon and lat internally - - ds = ds.rename({"gphif": "lat", "glamf": "lon"}) # TODO: Logging message about rename ds["grid"] = xr.DataArray( 0, attrs=sgrid.Grid2DMetadata( cf_role="grid_topology", topology_dimension=2, node_dimensions=("x", "y"), - node_coordinates=("lon", "lat"), + node_coordinates=("glamf", "gphif"), face_dimensions=( sgrid.DimDimPadding("x_center", "x", sgrid.Padding.LOW), sgrid.DimDimPadding("y_center", "y", sgrid.Padding.LOW), @@ -264,6 +261,9 @@ def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset | xr.Dat ) # NEMO models are always in degrees - ds["lon"].attrs["units"] = "degrees" - ds["lat"].attrs["units"] = "degrees" + ds["glamf"].attrs["units"] = "degrees" + ds["gphif"].attrs["units"] = "degrees" + + # Update to use lon and lat for internal naming + ds = sgrid.rename(ds, {"gphif": "lat", "glamf": "lon"}) # TODO: Logging message about rename return ds From 98395a0c12ef74897223705e77e99d29aa3ef730 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 13 Jan 2026 13:27:43 +0100 Subject: [PATCH 36/38] Remove setting of converters --- src/parcels/_core/fieldset.py | 6 ------ tests/test_advection.py | 7 +------ 2 files changed, 1 insertion(+), 12 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 9288ca8ff..9898463b5 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -10,7 +10,6 @@ import xarray as xr import xgcm -from parcels._core.converters import Geographic, GeographicPolar from parcels._core.field import Field, VectorField from parcels._core.utils import sgrid from parcels._core.utils.string import _assert_str_and_python_varname @@ -384,11 +383,6 @@ def from_sgrid_conventions( fields[varname] = Field(varname, ds[varname], grid, XLinear) fieldset = cls(list(fields.values())) - if mesh == "spherical": - if "U" in fieldset.fields: - fieldset.U.units = GeographicPolar() - if "V" in fieldset.fields: - fieldset.V.units = Geographic() return fieldset diff --git a/tests/test_advection.py b/tests/test_advection.py index dc5ae347c..4f7c64c58 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -438,12 +438,7 @@ def test_nemo_curvilinear_fieldset(): "kernel", [ AdvectionRK4, - pytest.param( - AdvectionRK4_3D, - marks=pytest.mark.xfail( - reason="from_nemo had 'fieldset.V.units = GeographicPolar()', I'm not sure _why_ this code is needed to get this to pass. To be further investigated.", - ), - ), + AdvectionRK4_3D, ], ) def test_nemo_3D_curvilinear_fieldset(kernel): From f726363febd7a4e53606f7243d5affeabe9cc345 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 13 Jan 2026 17:27:10 +0100 Subject: [PATCH 37/38] Select vector_interp_method based on U and V --- src/parcels/_core/fieldset.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 63e071238..26565929f 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -22,6 +22,7 @@ from parcels._typing import Mesh from parcels.convert import _discover_U_and_V, _ds_rename_using_standard_names, _maybe_rename_coords from parcels.interpolators import ( + CGrid_Velocity, Ux_Velocity, UxPiecewiseConstantFace, UxPiecewiseLinearNode, @@ -365,14 +366,15 @@ def from_sgrid_conventions( fields = {} if "U" in ds.data_vars and "V" in ds.data_vars: + vector_interp_method = XLinear_Velocity if _is_agrid(ds) else CGrid_Velocity fields["U"] = Field("U", ds["U"], grid, XLinear) fields["V"] = Field("V", ds["V"], grid, XLinear) - fields["UV"] = VectorField("UV", fields["U"], fields["V"], vector_interp_method=XLinear_Velocity) + fields["UV"] = VectorField("UV", fields["U"], fields["V"], vector_interp_method=vector_interp_method) if "W" in ds.data_vars: fields["W"] = Field("W", ds["W"], grid, XLinear) fields["UVW"] = VectorField( - "UVW", fields["U"], fields["V"], fields["W"], vector_interp_method=XLinear_Velocity + "UVW", fields["U"], fields["V"], fields["W"], vector_interp_method=vector_interp_method ) for varname in set(ds.data_vars) - set(fields.keys()) - skip_vars: @@ -548,6 +550,12 @@ def _get_mesh_type_from_sgrid_dataset(ds_sgrid: xr.Dataset) -> Mesh: return "spherical" if _is_coordinate_in_degrees(ds_sgrid[fpoint_x]) else "flat" +def _is_agrid(ds: xr.Dataset) -> bool: + # check if U and V are defined on the same dimensions + # if yes, interpret as A grid + return set(ds["U"].dims) == set(ds["V"].dims) + + def _is_coordinate_in_degrees(da: xr.DataArray) -> bool: match da.attrs.get("units"): case None: From 3e8fc759baafff96d633befbc4508bbb65264ace Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 13 Jan 2026 17:29:44 +0100 Subject: [PATCH 38/38] Update parcels.convert.nemo_to_sgrid args Remove **fields and replace with normal arg --- docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb | 2 +- src/parcels/convert.py | 2 +- tests/test_advection.py | 4 ++-- tests/test_convert.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index 007d50d1d..ba26f3558 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb @@ -73,7 +73,7 @@ "\n", "ds_coords = xr.open_dataset(data_folder / \"mesh_mask.nc4\", decode_times=False)\n", "ds_fset = parcels.convert.nemo_to_sgrid(\n", - " U=ds_fields[\"U\"], V=ds_fields[\"V\"], coords=ds_coords\n", + " fields=dict(U=ds_fields[\"U\"], V=ds_fields[\"V\"]), coords=ds_coords\n", ")\n", "\n", "fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)" diff --git a/src/parcels/convert.py b/src/parcels/convert.py index d2c7533dd..554ef8629 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -163,7 +163,7 @@ def _discover_U_and_V(ds: xr.Dataset, cf_standard_names_fallbacks) -> xr.Dataset return ds -def nemo_to_sgrid(*, coords: xr.Dataset, **fields: dict[str, xr.Dataset | xr.DataArray]): +def nemo_to_sgrid(*, fields: dict[str, xr.Dataset | xr.DataArray], coords: xr.Dataset): # TODO: Update docstring """Create a FieldSet from a xarray.Dataset from NEMO netcdf files. diff --git a/tests/test_advection.py b/tests/test_advection.py index fad137523..671c78d71 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -449,7 +449,7 @@ def test_nemo_curvilinear_fieldset(): V = xr.open_mfdataset(data_folder.glob("*V.nc4")) coords = xr.open_dataset(data_folder / "mesh_mask.nc4") - ds = parcels.convert.nemo_to_sgrid(U=U, V=V, coords=coords) + ds = parcels.convert.nemo_to_sgrid(fields=dict(U=U, V=V), coords=coords) fieldset = parcels.FieldSet.from_sgrid_conventions(ds) @@ -477,7 +477,7 @@ def test_nemo_3D_curvilinear_fieldset(kernel): W = xr.open_mfdataset(data_folder.glob("*W.nc")) coords = xr.open_dataset(data_folder / "coordinates.nc", decode_times=False) - ds = parcels.convert.nemo_to_sgrid(U=U["uo"], V=V["vo"], W=W["wo"], coords=coords) + ds = parcels.convert.nemo_to_sgrid(fields=dict(U=U["uo"], V=V["vo"], W=W["wo"]), coords=coords) fieldset = parcels.FieldSet.from_sgrid_conventions(ds) diff --git a/tests/test_convert.py b/tests/test_convert.py index 11a7029dc..eaf39d810 100644 --- a/tests/test_convert.py +++ b/tests/test_convert.py @@ -11,7 +11,7 @@ def test_nemo_to_sgrid(): V = xr.open_mfdataset(data_folder.glob("*V.nc4")) coords = xr.open_dataset(data_folder / "mesh_mask.nc4") - ds = convert.nemo_to_sgrid(U=U, V=V, coords=coords) + ds = convert.nemo_to_sgrid(fields=dict(U=U, V=V), coords=coords) assert ds["grid"].attrs == { "cf_role": "grid_topology",