Merge branch 'v4-dev' into support_1d_fields

Author: erikvansebille

.github/workflows/ci.yml

@@ -28,7 +28,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu, mac, windows]
+        os: [ubuntu, windows]
         pixi-environment: [test-latest]
         include:
           - os: ubuntu

.pre-commit-config.yaml

@@ -10,7 +10,7 @@ repos:
         types: [text]
         files: \.(json|ipynb)$
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.14.2
+    rev: v0.14.4
     hooks:
       - id: ruff
         name: ruff lint (.py)

CLAUDE.md

@@ -0,0 +1,10 @@
+## GitHub Interaction Guidelines
+
+- **NEVER impersonate the user on GitHub**, always sign off with something like
+  "[This is Claude Code on behalf of Jane Doe]"
+- Never create issues nor pull requests on the GitHub repository unless
+  explicitly instructed
+- Never post "update" messages, progress reports, or explanatory comments on
+  GitHub issues/PRs unless specifically instructed
+- When creating commits, always include a co-authorship trailer:
+  `Co-authored-by: Claude <[email protected]>`

docs/development/policies.md

@@ -1,5 +1,17 @@
 # Policies
 
+## Use of AI in development
+
+Many developers use AI Large Language Models to help them in their work. These LLMs have received both praise and criticism when it comes to software development.
+
+We accept that Parcels developers have their own motivation for using (or not using) AI. However, we have one policy that we expect all Parcels developers to follow:
+
+> It is ultimately your responsibility to understand the code that you commit.
+
+Remember that reviews are done by human maintainers - asking us to review code that an AI wrote but you don't understand isn't kind to these maintainers.
+
+The [CLAUDE.md](/CLAUDE.md) file in the repository has additional instructions for AI agents to follow when contributing to Parcels.
+
 ## Versioning
 
 Parcels follows [Intended Effort Versioning (EffVer)](https://jacobtomlinson.dev/effver/), where the version number (e.g., v2.1.0) is thought of as `MACRO.MESO.MICRO`.

src/parcels/_core/field.py

@@ -26,8 +26,6 @@
 from parcels._reprs import default_repr
 from parcels._typing import VectorType
 from parcels.interpolators import (
-    UXPiecewiseLinearNode,
-    XLinear,
     ZeroInterpolator,
     ZeroInterpolator_Vector,
 )
@@ -51,17 +49,11 @@ def _deal_with_errors(error, key, vector_type: VectorType):
         return 0
 
 
-_DEFAULT_INTERPOLATOR_MAPPING = {
-    XGrid: XLinear,
-    UxGrid: UXPiecewiseLinearNode,
-}
-
-
 class Field:
     """The Field class that holds scalar field data.
     The `Field` object is a wrapper around a xarray.DataArray or uxarray.UxDataArray object.
     Additionally, it holds a dynamic Callable procedure that is used to interpolate the field data.
-    During initialization, the user can supply a custom interpolation method that is used to interpolate the field data,
+    During initialization, the user is required to supply a custom interpolation method that is used to interpolate the field data,
     so long as the interpolation method has the correct signature.
 
     Notes
@@ -96,7 +88,7 @@ def __init__(
         name: str,
         data: xr.DataArray | ux.UxDataArray,
         grid: UxGrid | XGrid,
-        interp_method: Callable | None = None,
+        interp_method: Callable,
     ):
         if not isinstance(data, (ux.UxDataArray, xr.DataArray)):
             raise ValueError(
@@ -136,11 +128,8 @@ def __init__(
             raise e
 
         # Setting the interpolation method dynamically
-        if interp_method is None:
-            self._interp_method = _DEFAULT_INTERPOLATOR_MAPPING[type(self.grid)]
-        else:
-            assert_same_function_signature(interp_method, ref=ZeroInterpolator, context="Interpolation")
-            self._interp_method = interp_method
+        assert_same_function_signature(interp_method, ref=ZeroInterpolator, context="Interpolation")
+        self._interp_method = interp_method
 
         self.igrid = -1  # Default the grid index to -1
 

src/parcels/_core/fieldset.py

@@ -6,6 +6,7 @@
 
 import cf_xarray  # noqa: F401
 import numpy as np
+import uxarray as ux
 import xarray as xr
 import xgcm
 
@@ -14,10 +15,11 @@
 from parcels._core.utils.string import _assert_str_and_python_varname
 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
-from parcels.interpolators import XConstantField
+from parcels.interpolators import UXPiecewiseConstantFace, UXPiecewiseLinearNode, XConstantField, XLinear
 
 if TYPE_CHECKING:
     from parcels._core.basegrid import BaseGrid
@@ -238,22 +240,62 @@ def from_copernicusmarine(ds: xr.Dataset):
 
         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"], grid)
+            fields["U"] = Field("U", ds["U"], grid, XLinear)
+            fields["V"] = Field("V", ds["V"], grid, XLinear)
             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)
+                fields["W"] = Field("W", ds["W"], grid, XLinear)
                 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)
+            fields[varname] = Field(varname, ds[varname], grid, XLinear)
+
+        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.
+        """
+        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, _select_uxinterpolator(ds["U"]))
+            fields["V"] = Field("V", ds["V"], grid, _select_uxinterpolator(ds["U"]))
+            fields["U"].units = GeographicPolar()
+            fields["V"].units = Geographic()
+
+            if "W" in ds.data_vars:
+                fields["W"] = Field("W", ds["W"], grid, _select_uxinterpolator(ds["U"]))
+                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, _select_uxinterpolator(ds[varname]))
 
         return FieldSet(list(fields.values()))
 
@@ -365,11 +407,86 @@ 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:
+    # Common variable names for U and V found in UxDatasets
+    common_fesom_UV = [("unod", "vnod"), ("u", "v")]
+    common_fesom_W = ["w"]
+
+    if "W" not in ds:
+        for common_W in common_fesom_W:
+            if common_W in ds:
+                ds = _ds_rename_using_standard_names(ds, {common_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 common_U, common_V in common_fesom_UV:
+        if common_U in ds:
+            if common_V not in ds:
+                raise ValueError(
+                    f"Dataset has variable with standard name {common_U!r}, "
+                    f"but not the matching variable with standard name {common_V!r}. "
+                    "Please rename the appropriate variables in your dataset to have both 'U' and 'V' for Parcels simulation."
+                )
+            else:
+                ds = _ds_rename_using_standard_names(ds, {common_U: "U", common_V: "V"})
+                break
+
+        else:
+            if common_V in ds:
+                raise ValueError(
+                    f"Dataset has variable with standard name {common_V!r}, "
+                    f"but not the matching variable with standard name {common_U!r}. "
+                    "Please rename the appropriate variables in your dataset to have both 'U' and 'V' for Parcels simulation."
+                )
+            continue
+
+    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 = {
+        # (nz1,n_face): face-center laterally, layer centers vertically — piecewise constant
+        "nz1,n_face": UXPiecewiseConstantFace,
+        # (nz,n_node): node/corner laterally, layer interfaces vertically — barycentric lateral & linear vertical
+        "nz,n_node": UXPiecewiseLinearNode,
+    }
+    # Extract only spatial dimensions, neglecting time
+    da_spatial_dims = tuple(d for d in da.dims if d not in ("time",))
+    if len(da_spatial_dims) != 2:
+        raise ValueError(
+            "Fields on unstructured grids must have two spatial dimensions, one vertical (nz or nz1) and one lateral (n_face, n_edge, or n_node)"
+        )
+
+    # Construct key (string) for mapping to interpolator
+    # Find vertical and lateral tokens
+    vdim = None
+    ldim = None
+    for d in da_spatial_dims:
+        if d in ("nz", "nz1"):
+            vdim = d
+        if d in ("n_face", "n_node"):
+            ldim = d
+    # Map to supported interpolators
+    if vdim and ldim:
+        key = f"{vdim},{ldim}"
+        if key in supported_uxinterp_mapping.keys():
+            return supported_uxinterp_mapping[key]
+
+    return None

tests/test_advection.py

@@ -362,8 +362,8 @@ def test_stommelgyre_fieldset(kernel, rtol, grid_type):
     ds = stommel_gyre_dataset(grid_type=grid_type)
     grid = XGrid.from_dataset(ds)
     vector_interp_method = None if grid_type == "A" else CGrid_Velocity
-    U = Field("U", ds["U"], grid)
-    V = Field("V", ds["V"], grid)
+    U = Field("U", ds["U"], grid, interp_method=XLinear)
+    V = Field("V", ds["V"], grid, interp_method=XLinear)
     P = Field("P", ds["P"], grid, interp_method=XLinear)
     UV = VectorField("UV", U, V, vector_interp_method=vector_interp_method)
     fieldset = FieldSet([U, V, P, UV])
@@ -451,8 +451,8 @@ def test_nemo_curvilinear_fieldset():
     )
     grid = XGrid(xgcm_grid, mesh="spherical")
 
-    U = parcels.Field("U", ds["U"], grid)
-    V = parcels.Field("V", ds["V"], grid)
+    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 GoegraphicPolar for  C-Grid interpolation to work correctly
     UV = parcels.VectorField("UV", U, V, vector_interp_method=CGrid_Velocity)
@@ -536,9 +536,9 @@ def test_nemo_3D_curvilinear_fieldset(kernel):
     )
     grid = XGrid(xgcm_grid, mesh="spherical")
 
-    U = parcels.Field("U", ds["U"], grid)
-    V = parcels.Field("V", ds["V"], grid)
-    W = parcels.Field("W", ds["W"], grid)
+    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)

tests/test_field.py

@@ -9,37 +9,37 @@
 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 parcels.interpolators import UXPiecewiseConstantFace, UXPiecewiseLinearNode
+from parcels.interpolators import UXPiecewiseConstantFace, UXPiecewiseLinearNode, XLinear
 
 
 def test_field_init_param_types():
     data = datasets_structured["ds_2d_left"]
     grid = XGrid.from_dataset(data)
 
     with pytest.raises(TypeError, match="Expected a string for variable name, got int instead."):
-        Field(name=123, data=data["data_g"], grid=grid)
+        Field(name=123, data=data["data_g"], grid=grid, interp_method=XLinear)
 
     for name in ["a b", "123"]:
         with pytest.raises(
             ValueError,
             match=r"Received invalid Python variable name.*: not a valid identifier. HINT: avoid using spaces, special characters, and starting with a number.",
         ):
-            Field(name=name, data=data["data_g"], grid=grid)
+            Field(name=name, data=data["data_g"], grid=grid, interp_method=XLinear)
 
     with pytest.raises(
         ValueError,
         match=r"Received invalid Python variable name.*: it is a reserved keyword. HINT: avoid using the following names:.*",
     ):
-        Field(name="while", data=data["data_g"], grid=grid)
+        Field(name="while", data=data["data_g"], grid=grid, interp_method=XLinear)
 
     with pytest.raises(
         ValueError,
         match="Expected `data` to be a uxarray.UxDataArray or xarray.DataArray",
     ):
-        Field(name="test", data=123, grid=grid)
+        Field(name="test", data=123, grid=grid, interp_method=XLinear)
 
     with pytest.raises(ValueError, match="Expected `grid` to be a parcels UxGrid, or parcels XGrid"):
-        Field(name="test", data=data["data_g"], grid=123)
+        Field(name="test", data=data["data_g"], grid=123, interp_method=XLinear)
 
 
 @pytest.mark.parametrize(
@@ -66,6 +66,7 @@ def test_field_incompatible_combination(data, grid):
             name="test_field",
             data=data,
             grid=grid,
+            interp_method=XLinear,
         )
 
 
@@ -85,6 +86,7 @@ def test_field_init_structured_grid(data, grid):
         name="test_field",
         data=data,
         grid=grid,
+        interp_method=XLinear,
     )
     assert field.name == "test_field"
     assert field.data.equals(data)
@@ -113,6 +115,7 @@ def test_field_init_fail_on_float_time_dim():
             name="test_field",
             data=data,
             grid=grid,
+            interp_method=XLinear,
         )
 
 
@@ -128,7 +131,7 @@ def test_field_init_fail_on_float_time_dim():
 )
 def test_field_time_interval(data, grid):
     """Test creating a field."""
-    field = Field(name="test_field", data=data, grid=grid)
+    field = Field(name="test_field", data=data, grid=grid, interp_method=XLinear)
     assert field.time_interval.left == np.datetime64("2000-01-01")
     assert field.time_interval.right == np.datetime64("2001-01-01")
 
@@ -163,8 +166,8 @@ def invalid_interpolator_wrong_signature(particle_positions, grid_positions, inv
         return 0.0
 
     # Create component fields
-    U = Field(name="U", data=ds["data_g"], grid=grid)
-    V = Field(name="V", data=ds["data_g"], grid=grid)
+    U = Field(name="U", data=ds["data_g"], grid=grid, interp_method=XLinear)
+    V = Field(name="V", data=ds["data_g"], grid=grid, interp_method=XLinear)
 
     # Test invalid interpolator with wrong signature
     with pytest.raises(ValueError, match=".*incorrect name.*"):