diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 5864b9288..f0e48577f 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -289,6 +289,7 @@ def from_sgrid_conventions( ds: xr.Dataset, mesh: ptyping.TMesh | None = None, vector_fields: ptyping.VectorFields | NotSetType = NOTSET, + skip_field_data_validation: bool = False, ): # TODO: Update mesh to be discovered from the dataset metadata """Create a FieldSet from a dataset using SGRID convention metadata. @@ -307,6 +308,8 @@ def from_sgrid_conventions( Mapping of vector field names to tuples of component variable names in the dataset. For example, ``{"UV": ("U", "V"), "UVW": ("U", "V", "W")}``. If omitted (default), vector fields are auto-discovered from standard variable names (``U``/``V``/``W``). + skip_field_data_validation : bool, optional + If True, skip validation of the field data. This can be useful for performance reasons, but may lead to unexpected behavior if the field data is invalid. Default is False. Returns ------- @@ -321,7 +324,9 @@ def from_sgrid_conventions( See https://sgrid.github.io/sgrid/ for more information on the SGRID conventions. """ - model = StructuredModelData.from_sgrid_conventions(ds, mesh, vector_fields) + model = StructuredModelData.from_sgrid_conventions( + ds, mesh, vector_fields, skip_field_data_validation=skip_field_data_validation + ) return cls([model]) def describe(self, buf: IO | None = None) -> None: diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index 24542f7ba..9f31711d3 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -132,11 +132,19 @@ def preprocess_sgrid_model_data(ds: xr.Dataset) -> xr.Dataset: class StructuredModelData(ModelData): - def __init__(self, data: xr.Dataset, mesh: ptyping.TMesh, vector_field_components: ptyping.VectorFields): + def __init__( + self, + data: xr.Dataset, + mesh: ptyping.TMesh, + vector_field_components: ptyping.VectorFields, + skip_field_data_validation: bool = False, + ): if not isinstance(data, xr.Dataset): raise ValueError(f"Expected `data` to be an xarray.Dataset . Got {type(data)}") data = preprocess_sgrid_model_data(data) + if not skip_field_data_validation: + data = data.fillna(0) grid = XGrid(data, mesh) self.data = data @@ -181,7 +189,11 @@ def construct_fields(self) -> list[Field | VectorField]: @classmethod def from_sgrid_conventions( - cls, ds: xr.Dataset, mesh: ptyping.TMesh | None, vector_fields: ptyping.VectorFields | NotSetType + cls, + ds: xr.Dataset, + mesh: ptyping.TMesh | None, + vector_fields: ptyping.VectorFields | NotSetType, + skip_field_data_validation: bool = False, ) -> Self: ds = ds.copy() if mesh is None: @@ -212,7 +224,12 @@ def from_sgrid_conventions( vector_fields = resolve_vector_fields(ds, vector_fields) assert_valid_vector_fields(ds, vector_fields) - model = cls(ds, mesh=mesh, vector_field_components=vector_fields) + model = cls( + ds, + mesh=mesh, + vector_field_components=vector_fields, + skip_field_data_validation=skip_field_data_validation, + ) model._fields = model.construct_fields() for f in model._fields: if isinstance(f, Field): diff --git a/src/parcels/convert.py b/src/parcels/convert.py index c8cea9ce3..be2a756bb 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -541,6 +541,7 @@ def copernicusmarine_to_sgrid( ds.attrs.clear() # Clear global attributes from the merging 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"].data *= -1 diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 015d4110e..26310f703 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -330,6 +330,16 @@ def test_fieldset_from_sgrid_conventions(ds_name): assert len(fieldset.fields) > 0 +def test_fieldset_skip_field_data_validation(): + ds = datasets_structured["ds_2d_left"] + ds["U_A_grid"][:] = np.nan + + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") + assert np.isfinite(fieldset.U_A_grid.data).all() + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat", skip_field_data_validation=True) + assert np.isnan(fieldset.U_A_grid.data).all() + + def test_fieldset_add_error_on_duplicate_fields(): """Test that adding FieldSets with overlapping field names raises a ValueError.""" ds1 = datasets_structured["ds_2d_left"][["U_A_grid", "V_A_grid", "grid"]].rename({"U_A_grid": "U", "V_A_grid": "V"}) @@ -490,9 +500,9 @@ def test_fieldset_describe_backends(tmp_path): expected = """\ | Name | Type | Grid number | Interp method / value | Parcels backend | |:-------|:------------|--------------:|:------------------------|:------------------| -| U | Field | 0 | XLinear(...) | Zarr | -| V | Field | 0 | XLinear(...) | Zarr | -| W | Field | 0 | XLinear(...) | Zarr | +| U | Field | 0 | XLinear(...) | NumPy | +| V | Field | 0 | XLinear(...) | NumPy | +| W | Field | 0 | XLinear(...) | NumPy | | UV | VectorField | 0 | CGrid_Velocity(...) | - | | UVW | VectorField | 0 | CGrid_Velocity(...) | - |