From 0b20ab11d6448d1da5111c52a3e039db87f6253d Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 8 Oct 2025 16:25:12 -0400 Subject: [PATCH 1/5] Add method to import from fesom2 Add single test for stommel gyre generic dataset --- src/parcels/_core/fieldset.py | 89 ++++++++++++++++++++++++++++++++++- tests/test_fieldset.py | 9 ++++ 2 files changed, 97 insertions(+), 1 deletion(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 7cbe02499..d638cc289 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -6,13 +6,16 @@ import cf_xarray # noqa: F401 import numpy as np +import uxarray as ux import xarray as xr import xgcm +from parcels import interpolators from parcels._core.converters import Geographic, GeographicPolar from parcels._core.field import Field, VectorField from parcels._core.utils.time import get_datetime_type_calendar from parcels._core.utils.time import is_compatible as datetime_is_compatible +from parcels._core.uxgrid import UxGrid from parcels._core.xgrid import _DEFAULT_XGCM_KWARGS, XGrid from parcels._logger import logger from parcels._typing import Mesh @@ -251,6 +254,55 @@ def from_copernicusmarine(ds: xr.Dataset): return FieldSet(list(fields.values())) + def from_fesom2(ds: ux.UxDataset): + """Create a FieldSet from a FESOM2 uxarray.UxDataset. + + Parameters + ---------- + ds : uxarray.UxDataset + uxarray.UxDataset as obtained from the uxarray package. + + Returns + ------- + FieldSet + FieldSet object containing the fields from the dataset that can be used for a Parcels simulation. + + """ + location_to_interpmethod = { + "face": interpolators.UXPiecewiseConstantFace, + "node": interpolators.UXPiecewiseLinearNode, + } # TODO : Edge registered interpolator, independent vertical grid detection from lateral node placement + + ds = ds.copy() + ds_dims = list(ds.dims) + if not all(dim in ds_dims for dim in ["time", "nz", "nz1"]): + raise ValueError( + f"Dataset missing one of the required dimensions 'time', 'nz', or 'nz1'. Found dimensions {ds_dims}" + ) + grid = UxGrid(ds.uxgrid, z=ds.coords["nz"]) + ds = _discover_fesom2_U_and_V(ds) + + fields = {} + if "U" in ds.data_vars and "V" in ds.data_vars: + fields["U"] = Field("U", ds["U"], grid, location_to_interpmethod[ds["U"].attrs.get("location")]) + fields["V"] = Field("V", ds["V"], grid, location_to_interpmethod[ds["V"].attrs.get("location")]) + fields["U"].units = GeographicPolar() + fields["V"].units = Geographic() + + if "W" in ds.data_vars: + ds["W"] -= ds[ + "W" + ] # Negate W to convert from up positive to down positive (as that's the direction of positive z) + fields["W"] = Field("W", ds["W"], grid, location_to_interpmethod[ds["W"].attrs.get("location")]) + fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"]) + else: + fields["UV"] = VectorField("UV", fields["U"], fields["V"]) + + for varname in set(ds.data_vars) - set(fields.keys()): + fields[varname] = Field(varname, ds[varname], grid) + + return FieldSet(list(fields.values())) + class CalendarError(Exception): # TODO: Move to a parcels errors module """Exception raised when the calendar of a field is not compatible with the rest of the Fields. The user should ensure that they only add fields to a FieldSet that have compatible CFtime calendars.""" @@ -359,7 +411,42 @@ def _discover_copernicusmarine_U_and_V(ds: xr.Dataset) -> xr.Dataset: return ds -def _ds_rename_using_standard_names(ds: xr.Dataset, name_dict: dict[str, str]) -> xr.Dataset: +def _discover_fesom2_U_and_V(ds: ux.UxDataset) -> ux.UxDataset: + # Assumes that the dataset has U and V data + + cf_UV_standard_name_fallbacks = [("unod", "vnod"), ("u", "v")] + cf_W_standard_name_fallbacks = ["w"] + + if "W" not in ds: + for cf_standard_name_W in cf_W_standard_name_fallbacks: + 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_UV_standard_name_fallbacks: + 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 _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}) diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 85d082d69..d25d2fd62 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -11,6 +11,7 @@ from parcels._datasets.structured.circulation_models import datasets as datasets_circulation_models from parcels._datasets.structured.generic import T as T_structured from parcels._datasets.structured.generic import datasets as datasets_structured +from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from tests import utils ds = datasets_structured["ds_2d_left"] @@ -270,3 +271,11 @@ def test_fieldset_from_copernicusmarine_with_W(caplog): assert "UV" not in fieldset.fields assert "UVW" in fieldset.fields assert "renamed it to 'W'" in caplog.text + + +def test_fieldset_from_fesom2(): + ds = datasets_unstructured["stommel_gyre_delaunay"] + fieldset = FieldSet.from_fesom2(ds) + assert "U" in fieldset.fields + assert "V" in fieldset.fields + assert "UVW" in fieldset.fields From aad13770bf9ec4c65c1bd6dcebe46cf785469086 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 8 Oct 2025 17:05:51 -0400 Subject: [PATCH 2/5] Remove assignment of interpolator in from_fesom2 We don't want to provide this convenience here, due to the fact that some variable placements are not supported (e.g. on an edge), but we don't want to prevent someone who is building a new interpolator from leveraging this new from_fesom2 method. --- src/parcels/_core/fieldset.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index d638cc289..5eea4c448 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -282,6 +282,16 @@ def from_fesom2(ds: ux.UxDataset): grid = UxGrid(ds.uxgrid, z=ds.coords["nz"]) ds = _discover_fesom2_U_and_V(ds) + for var in ds.data_vars: + if "location" not in ds[var].attrs: + raise ValueError( + f"Variable {var!r} missing required attribute 'location' to indicate its placement on the mesh." + ) + if ds[var].attrs["location"] not in location_to_interpmethod: + raise ValueError( + f"Variable {var!r} has attribute 'location' with unexpected value {ds[var].attrs['location']!r}. Expected one of {list(location_to_interpmethod.keys())}" + ) + fields = {} if "U" in ds.data_vars and "V" in ds.data_vars: fields["U"] = Field("U", ds["U"], grid, location_to_interpmethod[ds["U"].attrs.get("location")]) From 2f5aa563f9b2a068d65dd1a6b9913bcf7becc5ab Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 8 Oct 2025 17:10:00 -0400 Subject: [PATCH 3/5] Remove assignment of interpolator in from_fesom2 We don't want to provide this convenience here, due to the fact that some variable placements are not supported (e.g. on an edge), but we don't want to prevent someone who is building a new interpolator from leveraging this new from_fesom2 method. --- src/parcels/_core/fieldset.py | 22 +++------------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 5eea4c448..b15dca27b 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 interpolators from parcels._core.converters import Geographic, GeographicPolar from parcels._core.field import Field, VectorField from parcels._core.utils.time import get_datetime_type_calendar @@ -268,11 +267,6 @@ def from_fesom2(ds: ux.UxDataset): FieldSet object containing the fields from the dataset that can be used for a Parcels simulation. """ - location_to_interpmethod = { - "face": interpolators.UXPiecewiseConstantFace, - "node": interpolators.UXPiecewiseLinearNode, - } # TODO : Edge registered interpolator, independent vertical grid detection from lateral node placement - ds = ds.copy() ds_dims = list(ds.dims) if not all(dim in ds_dims for dim in ["time", "nz", "nz1"]): @@ -282,20 +276,10 @@ def from_fesom2(ds: ux.UxDataset): grid = UxGrid(ds.uxgrid, z=ds.coords["nz"]) ds = _discover_fesom2_U_and_V(ds) - for var in ds.data_vars: - if "location" not in ds[var].attrs: - raise ValueError( - f"Variable {var!r} missing required attribute 'location' to indicate its placement on the mesh." - ) - if ds[var].attrs["location"] not in location_to_interpmethod: - raise ValueError( - f"Variable {var!r} has attribute 'location' with unexpected value {ds[var].attrs['location']!r}. Expected one of {list(location_to_interpmethod.keys())}" - ) - fields = {} if "U" in ds.data_vars and "V" in ds.data_vars: - fields["U"] = Field("U", ds["U"], grid, location_to_interpmethod[ds["U"].attrs.get("location")]) - fields["V"] = Field("V", ds["V"], grid, location_to_interpmethod[ds["V"].attrs.get("location")]) + fields["U"] = Field("U", ds["U"], grid) + fields["V"] = Field("V", ds["V"]) fields["U"].units = GeographicPolar() fields["V"].units = Geographic() @@ -303,7 +287,7 @@ def from_fesom2(ds: ux.UxDataset): ds["W"] -= ds[ "W" ] # Negate W to convert from up positive to down positive (as that's the direction of positive z) - fields["W"] = Field("W", ds["W"], grid, location_to_interpmethod[ds["W"].attrs.get("location")]) + fields["W"] = Field("W", ds["W"], grid) fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"]) else: fields["UV"] = VectorField("UV", fields["U"], fields["V"]) From 582656d8e4430c8ca993ad1f00ea78cfeecaef64 Mon Sep 17 00:00:00 2001 From: Joe Schoonover <11430768+fluidnumerics-joe@users.noreply.github.com> Date: Thu, 9 Oct 2025 06:19:50 -0400 Subject: [PATCH 4/5] Add grid parameters to V field creation Co-authored-by: Erik van Sebille --- src/parcels/_core/fieldset.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index b15dca27b..11291e053 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -279,7 +279,7 @@ def from_fesom2(ds: ux.UxDataset): fields = {} if "U" in ds.data_vars and "V" in ds.data_vars: fields["U"] = Field("U", ds["U"], grid) - fields["V"] = Field("V", ds["V"]) + fields["V"] = Field("V", ds["V"], grid) fields["U"].units = GeographicPolar() fields["V"].units = Geographic() From 25fa2d39d0859712d70a14a30280ba23c91827ec Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Thu, 9 Oct 2025 06:31:44 -0400 Subject: [PATCH 5/5] Remove sign change on vertical velocity --- src/parcels/_core/fieldset.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 11291e053..16e2772bb 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -284,9 +284,6 @@ def from_fesom2(ds: ux.UxDataset): fields["V"].units = Geographic() if "W" in ds.data_vars: - ds["W"] -= ds[ - "W" - ] # Negate W to convert from up positive to down positive (as that's the direction of positive z) fields["W"] = Field("W", ds["W"], grid) fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"]) else: