use zero for land-values throughout Parcels (#2762)

* Add _assert_no_nan_in_field_data function

* Use fillna in convert functions

* Adding fill_nan_to_zero and assert_valid_fields options

* Changing to one skip_field_data_validation option

Author: erikvansebille

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:

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):

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

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(...)     | -                 |