Merge branch 'v4-dev' into MyST-quickstart

Author: reint-fischer

.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
 
@@ -217,6 +206,9 @@ def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True):
             _ei = None
         else:
             _ei = particles.ei[:, self.igrid]
+        z = np.atleast_1d(z)
+        y = np.atleast_1d(y)
+        x = np.atleast_1d(x)
 
         particle_positions, grid_positions = _get_positions(self, time, z, y, x, particles, _ei)
 
@@ -300,6 +292,9 @@ def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True):
             _ei = None
         else:
             _ei = particles.ei[:, self.igrid]
+        z = np.atleast_1d(z)
+        y = np.atleast_1d(y)
+        x = np.atleast_1d(x)
 
         particle_positions, grid_positions = _get_positions(self.U, time, z, y, x, particles, _ei)
 

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,9 +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 UXPiecewiseConstantFace, UXPiecewiseLinearNode, XConstantField, XLinear
 
 if TYPE_CHECKING:
     from parcels._core.basegrid import BaseGrid
@@ -116,7 +119,7 @@ def add_field(self, field: Field, name: str | None = None):
 
         self.fields[name] = field
 
-    def add_constant_field(self, name: str, value, mesh: Mesh = "flat"):
+    def add_constant_field(self, name: str, value, mesh: Mesh = "spherical"):
         """Wrapper function to add a Field that is constant in space,
            useful e.g. when using constant horizontal diffusivity
 
@@ -134,16 +137,15 @@ def add_constant_field(self, name: str, value, mesh: Mesh = "flat"):
                correction for zonal velocity U near the poles.
             2. flat: No conversion, lat/lon are assumed to be in m.
         """
-        ds = xr.Dataset({name: (["time", "lat", "lon", "depth"], np.full((1, 1, 1, 1), value))})
-        grid = XGrid(xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS))
-        self.add_field(
-            Field(
-                name,
-                ds[name],
-                grid,
-                interp_method=None,  # TODO : Need to define an interpolation method for constants
-            )
+        ds = xr.Dataset(
+            {name: (["lat", "lon"], np.full((1, 1), value))},
+            coords={"lat": (["lat"], [0], {"axis": "Y"}), "lon": (["lon"], [0], {"axis": "X"})},
+        )
+        xgrid = xgcm.Grid(
+            ds, coords={"X": {"left": "lon"}, "Y": {"left": "lat"}}, autoparse_metadata=False, **_DEFAULT_XGCM_KWARGS
         )
+        grid = XGrid(xgrid, mesh=mesh)
+        self.add_field(Field(name, ds[name], grid, interp_method=XConstantField))
 
     def add_constant(self, name, value):
         """Add a constant to the FieldSet. Note that all constants are
@@ -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, 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)
+            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

src/parcels/_core/xgrid.py

@@ -289,22 +289,13 @@ def search(self, z, y, x, ei=None):
         else:
             zi, zeta = np.zeros(z.shape, dtype=int), np.zeros(z.shape, dtype=float)
 
-        if ds.lon.ndim == 1:
-            yi, eta = _search_1d_array(ds.lat.values, y)
-            xi, xsi = _search_1d_array(ds.lon.values, x)
-            return {
-                "Z": {"index": zi, "bcoord": zeta},
-                "Y": {"index": yi, "bcoord": eta},
-                "X": {"index": xi, "bcoord": xsi},
-            }
+        if "X" in self.axes and "Y" in self.axes and ds.lon.ndim == 2:
+            yi, xi = None, None
+            if ei is not None:
+                axis_indices = self.unravel_index(ei)
+                xi = axis_indices.get("X")
+                yi = axis_indices.get("Y")
 
-        yi, xi = None, None
-        if ei is not None:
-            axis_indices = self.unravel_index(ei)
-            xi = axis_indices.get("X")
-            yi = axis_indices.get("Y")
-
-        if ds.lon.ndim == 2:
             yi, eta, xi, xsi = _search_indices_curvilinear_2d(self, y, x, yi, xi)
 
             return {
@@ -313,7 +304,24 @@ def search(self, z, y, x, ei=None):
                 "X": {"index": xi, "bcoord": xsi},
             }
 
-        raise NotImplementedError("Searching in >2D lon/lat arrays is not implemented yet.")
+        if "X" in self.axes and ds.lon.ndim > 2:
+            raise NotImplementedError("Searching in >2D lon/lat arrays is not implemented yet.")
+
+        if "Y" in self.axes:
+            yi, eta = _search_1d_array(ds.lat.values, y)
+        else:
+            yi, eta = np.zeros(y.shape, dtype=int), np.zeros(y.shape, dtype=float)
+
+        if "X" in self.axes:
+            xi, xsi = _search_1d_array(ds.lon.values, x)
+        else:
+            xi, xsi = np.zeros(x.shape, dtype=int), np.zeros(x.shape, dtype=float)
+
+        return {
+            "Z": {"index": zi, "bcoord": zeta},
+            "Y": {"index": yi, "bcoord": eta},
+            "X": {"index": xi, "bcoord": xsi},
+        }
 
     @cached_property
     def _fpoint_info(self):