From d7bc0f229eba6afcf78d473963a1497f93cfeb0a Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 4 Aug 2026 17:02:42 +0800 Subject: [PATCH 01/13] Add design spec for removing xgcm dependency --- .../specs/2026-08-04-remove-xgcm-design.md | 129 ++++++++++++++++++ 1 file changed, 129 insertions(+) create mode 100644 docs/superpowers/specs/2026-08-04-remove-xgcm-design.md diff --git a/docs/superpowers/specs/2026-08-04-remove-xgcm-design.md b/docs/superpowers/specs/2026-08-04-remove-xgcm-design.md new file mode 100644 index 000000000..310b53961 --- /dev/null +++ b/docs/superpowers/specs/2026-08-04-remove-xgcm-design.md @@ -0,0 +1,129 @@ +# Remove xgcm Dependency — Design Spec + +**Date:** 2026-08-04 +**Branch:** remove-xgcm +**Approach:** Option B — full SGRID-native refactor (remove XgcmLike* adapter layer) + +## Motivation + +All grid topology information previously sourced from xgcm (axis directions, staggering positions, face/node relationships) is now fully expressed in SGRID metadata attached to the dataset. The XgcmLike* adapter layer was always a temporary bridge. Removing xgcm simplifies the dependency graph, removes COMODO metadata coupling, and lets the codebase work entirely with the standardised SGRID model. + +--- + +## Section 1: New Position Vocabulary + +Replace xgcm position strings (`"center"`, `"left"`, `"right"`, `"inner"`, `"outer"`) with a new `GridPosition` type in `_typing.py`: + +```python +from parcels._sgrid.core import Padding + +GridPosition = Literal["face"] | Padding +``` + +| Old (xgcm) | New (SGRID) | +| ---------- | -------------- | +| `"center"` | `"face"` | +| `"right"` | `Padding.LOW` | +| `"left"` | `Padding.HIGH` | +| `"inner"` | `Padding.BOTH` | +| `"outer"` | `Padding.NONE` | + +**Removed from `_typing.py`:** + +- `XgcmAxisPosition` +- `XgcmAxes` +- `import xgcm` (TYPE_CHECKING block) + +**Removed from `_sgrid/core.py`:** + +- `SGRID_PADDING_TO_XGCM_POSITION` dict +- `xgcm_parse_sgrid()` function + +**Added to `_sgrid/accessor.py`:** + +```python +def get_dim_position(grid: SGrid2DMetadata, dim: str) -> GridPosition: + """Returns 'face' or the Padding value for a given dimension.""" +``` + +This replaces `get_xgcm_position_from_dim_name`. Uses the existing `_get_axis_info` helper internally. + +--- + +## Section 2: `xgrid.py` Changes + +### Removed entirely + +- `XgcmLikeAxis` dataclass +- `XgcmLikeGrid` class +- `construct_xgcm_axes_object` function +- `self.xgcm_grid` attribute on `XGrid` +- `_DEFAULT_XGCM_KWARGS` +- `import xgcm` and `import xgcm.axis` (TYPE_CHECKING blocks) +- Import of `SGRID_PADDING_TO_XGCM_POSITION` + +`self.sgrid_metadata` (already stored on `XGrid`) becomes the sole source of grid topology truth. + +### Changed function signatures + +| Function | Old signature | New signature | +| -------------------------------------------- | ----------------------------------------------------------------------------- | ------------------------------------------------------- | +| `get_cell_count_along_dim` | `(ds, axis: xgcm.axis.Axis)` | `(ds, fnp: FaceNodePadding)` — uses `ds[fnp.face].size` | +| `get_time` | `(ds, axis: xgcm.axis.Axis)` | `(ds, time_dim: str)` — uses `ds[time_dim].values` | +| `_get_xgrid_axes` | `(grid: xgcm.Grid)` | `(metadata: SGrid2DMetadata)` | +| `assert_all_field_dims_have_axis` | `(da, xgcm_grid: xgcm.Grid)` | `(da, metadata: SGrid2DMetadata)` | +| `assert_valid_lat_lon` | `(da_lat, da_lon, axes: XgcmAxes)` | `(da_lat, da_lon, metadata: SGrid2DMetadata)` | +| `assert_all_dimensions_correspond_with_axis` | `(da, axes: XgcmAxes)` | `(da, metadata: SGrid2DMetadata)` | +| `assert_valid_field_array` | `(da, axes: XgcmAxes)` | `(da, metadata: SGrid2DMetadata)` | +| `_convert_center_pos_to_fpoint` | `(xgcm_position: XgcmAxisPosition, f_points_xgcm_position: XgcmAxisPosition)` | `(position: GridPosition, f_point_position: Padding)` | + +### Internal logic changes + +- `get_axis_from_dim_name(axes, dim)` → replaced by `_get_dim_to_axis_mapping(metadata).get(dim)` from accessor +- `get_xgcm_position_from_dim_name(axes, dim)` → replaced by `get_dim_position(metadata, dim)` from accessor +- `_fpoint_info` → returns `dict[XgridAxis, Padding]` (was `dict[XgridAxis, str]` with xgcm strings) +- `localize()` → uses `self.sgrid_metadata` instead of `self.xgcm_grid` +- `get_axis_dim_mapping()` → uses `_get_dim_to_axis_mapping(self.sgrid_metadata)` directly +- `_convert_center_pos_to_fpoint`: `"center"` branch becomes `"face"` check; `"inner"/"right"` checks become `Padding.BOTH / Padding.LOW` +- `XGrid.lon/lat/depth/_datetimes` — currently check `self.xgcm_grid.axes["X"/"Y"/"Z"/"T"]` to determine axis presence; replaced by checking `_get_dim_to_axis_mapping(self.sgrid_metadata)` for spatial axes, and `"time" in self._ds.dims` for the time axis + +--- + +## Section 3: Test Changes + +### `tests/datasets/test_structured.py` + +Replace `xgcm.Grid` calls with SGRID-native checks: + +```python +# Old +grid = xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS) +for _axis_name, axis in grid.axes.items(): + for pos, _dim_name in axis.coords.items(): + assert pos in ["left", "center"] + +# New +metadata = ds.sgrid.metadata +for fnp in metadata.face_dimensions: + assert get_dim_position(metadata, fnp.face) == "face" + assert get_dim_position(metadata, fnp.node) in (Padding.HIGH, Padding.LOW, Padding.BOTH, Padding.NONE) +``` + +The specific padding assertion per test depends on the dataset fixture (`ds_2d_left` vs `ds_2d_right`). + +### `tests/sgrid/test_sgrid.py` + +- Delete the two tests at lines 259–285 that call `xgcm_parse_sgrid()` and `xgcm.Grid(...)` — they were testing the now-removed bridge function +- Remove `import xgcm` and `SGRID_PADDING_TO_XGCM_POSITION` from imports + +### `pyproject.toml` + +Remove `"xgcm >=0.9.0"` from the dependencies list. + +--- + +## Out of Scope + +- Changing interpolation logic or search algorithms in `xgrid.py` +- Renaming `XgcmAxisDirection` / `CfAxis` type aliases (they don't reference xgcm at runtime) +- Any changes to SGRID parsing or accessor logic beyond adding `get_dim_position` From 218f5a5c3822ef10cbf839309a2d7fb0d61f0a9d Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 4 Aug 2026 17:15:21 +0800 Subject: [PATCH 02/13] Add implementation plan for removing xgcm dependency --- .../plans/2026-08-04-remove-xgcm.md | 896 ++++++++++++++++++ 1 file changed, 896 insertions(+) create mode 100644 docs/superpowers/plans/2026-08-04-remove-xgcm.md diff --git a/docs/superpowers/plans/2026-08-04-remove-xgcm.md b/docs/superpowers/plans/2026-08-04-remove-xgcm.md new file mode 100644 index 000000000..0d52f463b --- /dev/null +++ b/docs/superpowers/plans/2026-08-04-remove-xgcm.md @@ -0,0 +1,896 @@ +# Remove xgcm Dependency Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** Remove xgcm as a runtime dependency by replacing the XgcmLike adapter layer with direct use of SGRID metadata throughout `xgrid.py`. + +**Architecture:** Three sequential tasks: (1) add new SGRID-native types and helpers as pure additions, (2) refactor xgrid.py to use them instead of the adapter layer, (3) delete the now-dead bridge functions and fix the two test files that called xgcm directly. + +**Tech Stack:** Python 3.11+, xarray, parcels SGRID metadata (`SGrid2DMetadata`, `FaceNodePadding`, `Padding`) + +## Global Constraints + +- Never import `xgcm` outside of `TYPE_CHECKING` blocks — and those blocks are being removed too +- Preserve all existing function names (they are imported by `model.py` and may be used elsewhere) +- All commits must include `Co-authored-by: Claude ` +- Run `pytest tests/sgrid/ tests/datasets/ -x` after each task before committing + +--- + +## File Map + +| File | Change | +| ----------------------------------- | -------------------------------------------------------------------------------------- | +| `src/parcels/_typing.py` | Add `GridPosition`; remove `XgcmAxisPosition`, `XgcmAxes`, xgcm TYPE_CHECKING import | +| `src/parcels/_sgrid/accessor.py` | Add `get_dim_position()` | +| `src/parcels/_core/xgrid.py` | Major refactor — remove adapter classes, update all functions to use `SGrid2DMetadata` | +| `src/parcels/_sgrid/core.py` | Remove `SGRID_PADDING_TO_XGCM_POSITION` and `xgcm_parse_sgrid()` | +| `src/parcels/_sgrid/__init__.py` | Remove `xgcm_parse_sgrid` from imports and `__all__` | +| `tests/sgrid/test_sgrid.py` | Remove xgcm import + two xgcm tests; add `test_get_dim_position` | +| `tests/datasets/test_structured.py` | Replace `xgcm.Grid` calls with SGRID-native assertions | +| `pyproject.toml` | Remove `"xgcm >=0.9.0"` | + +--- + +### Task 1: Add `GridPosition` type and `get_dim_position()` helper + +Pure additions — nothing is deleted. All existing code continues to work. + +**Files:** + +- Modify: `src/parcels/_typing.py` +- Modify: `src/parcels/_sgrid/accessor.py` +- Modify: `tests/sgrid/test_sgrid.py` + +**Interfaces:** + +- Produces: `GridPosition = Literal["face"] | Padding` in `_typing.py` +- Produces: `get_dim_position(grid: SGrid2DMetadata, dim: str) -> Literal["face"] | Padding` in `_sgrid/accessor.py` + +- [ ] **Step 1: Write failing tests for `get_dim_position`** + +Add to `tests/sgrid/test_sgrid.py` (after the existing imports, before existing test functions): + +```python +from parcels._sgrid.accessor import get_dim_position + + +def test_get_dim_position_face_dims(): + """Face dimensions return 'face'.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) + # face_dimensions = (FaceNodePadding("face_dimension1", "node_dimension1", Padding.LOW), ...) + assert get_dim_position(metadata, "face_dimension1") == "face" + assert get_dim_position(metadata, "face_dimension2") == "face" + + +def test_get_dim_position_node_dims(): + """Node dimensions return their Padding value.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) + assert get_dim_position(metadata, "node_dimension1") == sgrid.Padding.LOW + assert get_dim_position(metadata, "node_dimension2") == sgrid.Padding.LOW + + +def test_get_dim_position_vertical(): + """Vertical face and node dimensions are handled.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=True, with_node_coordinates=False) + # vertical_dimensions = (FaceNodePadding("vertical_dimensions_dim1", "vertical_dimensions_dim2", Padding.LOW),) + assert get_dim_position(metadata, "vertical_dimensions_dim1") == "face" + assert get_dim_position(metadata, "vertical_dimensions_dim2") == sgrid.Padding.LOW + + +def test_get_dim_position_unknown_dim_raises(): + """Unknown dimensions raise ValueError.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) + with pytest.raises(ValueError, match="not a spatial SGRID dimension"): + get_dim_position(metadata, "nonexistent_dim") +``` + +- [ ] **Step 2: Run tests to confirm they fail** + +```bash +pytest tests/sgrid/test_sgrid.py::test_get_dim_position_face_dims -x -v +``` + +Expected: `ImportError` or `AttributeError` — `get_dim_position` does not exist yet. + +- [ ] **Step 3: Add `GridPosition` to `_typing.py`** + +In `src/parcels/_typing.py`, after the existing imports block, add: + +```python +from parcels._sgrid.core import Padding +``` + +Then after `XgcmAxisDirection = CfAxisSpatial | Literal["T"]`, add: + +```python +GridPosition = Literal["face"] | Padding +``` + +No existing lines are removed yet — that happens in Task 2. + +- [ ] **Step 4: Add `get_dim_position` to `accessor.py`** + +In `src/parcels/_sgrid/accessor.py`, add this function after `_get_axis_info`: + +```python +def get_dim_position(grid: SGrid2DMetadata, dim: str) -> "Literal['face'] | Padding": + """Returns 'face' if dim is a face dimension, or the Padding value if it is a node dimension. + + Replaces xgcm's position string vocabulary ('center', 'left', 'right', 'inner', 'outer') + with SGRID-native types. + """ + axis_info = _get_axis_info(grid) + if dim not in axis_info: + raise ValueError(f"Dimension {dim!r} is not a spatial SGRID dimension in this grid.") + fnp, is_node = axis_info[dim] + return fnp.padding if is_node else "face" +``` + +Add the `Literal` import to the existing `from typing import Any, Literal, cast` line in accessor.py (it already has `Literal` — confirm and leave as-is if so). + +- [ ] **Step 5: Run tests to confirm they pass** + +```bash +pytest tests/sgrid/test_sgrid.py::test_get_dim_position_face_dims tests/sgrid/test_sgrid.py::test_get_dim_position_node_dims tests/sgrid/test_sgrid.py::test_get_dim_position_vertical tests/sgrid/test_sgrid.py::test_get_dim_position_unknown_dim_raises -v +``` + +Expected: 4 PASSED. + +- [ ] **Step 6: Run full test subset to confirm no regressions** + +```bash +pytest tests/sgrid/ tests/datasets/ -x -q +``` + +Expected: all pass (same as before this task). + +- [ ] **Step 7: Commit** + +```bash +git add src/parcels/_typing.py src/parcels/_sgrid/accessor.py tests/sgrid/test_sgrid.py +git commit -m "feat: add GridPosition type and get_dim_position() SGRID helper + +Co-authored-by: Claude " +``` + +--- + +### Task 2: Refactor `xgrid.py` to use SGRID metadata directly + +Remove `XgcmLikeAxis`, `XgcmLikeGrid`, `construct_xgcm_axes_object`, and `self.xgcm_grid`. Update every function in `xgrid.py` that referenced the adapter layer to use `SGrid2DMetadata` / `FaceNodePadding` / `Padding` directly. Also clean up the now-dead type aliases from `_typing.py`. + +**Files:** + +- Modify: `src/parcels/_core/xgrid.py` +- Modify: `src/parcels/_typing.py` + +**Interfaces:** + +- Consumes: `get_dim_position(grid, dim)` from Task 1 +- Consumes: `GridPosition` from Task 1 +- Consumes: `_get_dim_to_axis_mapping(metadata)` from `_sgrid/accessor.py` (already imported) +- Produces: All public function signatures updated to take `SGrid2DMetadata` instead of `XgcmAxes`/`xgcm.Grid` + +- [ ] **Step 1: Confirm baseline tests pass before touching anything** + +```bash +pytest tests/sgrid/ tests/datasets/ -x -q +``` + +Expected: all pass. If not, do not proceed — fix the cause first. + +- [ ] **Step 2: Rewrite `xgrid.py` imports and module-level constants** + +At the top of `src/parcels/_core/xgrid.py`, make these changes: + +Remove the `TYPE_CHECKING` block entirely: + +```python +# DELETE these lines: +if TYPE_CHECKING: + import xgcm.axis +``` + +Remove the `_DEFAULT_XGCM_KWARGS` constant: + +```python +# DELETE this line: +_DEFAULT_XGCM_KWARGS: dict[str, Any] = {"padding": "fill"} +``` + +Remove the import of `SGRID_PADDING_TO_XGCM_POSITION`: + +```python +# DELETE this line: +from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION +``` + +Add `get_dim_position` to the existing accessor import: + +```python +# Change this: +from parcels._sgrid.accessor import _get_dim_to_axis_mapping +# To: +from parcels._sgrid.accessor import _get_dim_to_axis_mapping, get_dim_position +``` + +Add `FaceNodePadding` and `Padding` to the sgrid imports. The existing `import parcels._sgrid as sgrid` is kept — use `sgrid.Padding` and `sgrid.FaceNodePadding` in the code below. + +Also update the `_typing` import to add `GridPosition`: + +```python +# The existing line: +import parcels._typing as ptyping +# No change to the line itself — access GridPosition as ptyping.GridPosition +``` + +Remove `Any` from the typing imports if it's only used by `_DEFAULT_XGCM_KWARGS` (check first — if used elsewhere, keep it). + +- [ ] **Step 3: Replace standalone helper functions** + +Replace `get_cell_count_along_dim` (lines 29-33): + +```python +def get_cell_count_along_dim(ds: xr.Dataset, fnp: sgrid.FaceNodePadding) -> int: + return ds[fnp.face].size - 1 +``` + +Replace `get_time` (lines 36-37): + +```python +def get_time(ds: xr.Dataset, time_dim: str) -> npt.NDArray: + return ds[time_dim].values +``` + +Replace `_get_xgrid_axes` (lines 40-42): + +```python +def _get_xgrid_axes(metadata: sgrid.SGrid2DMetadata, ds_dims: set[str]) -> list[ptyping.XgridAxis]: + dim_to_axis = _get_dim_to_axis_mapping(metadata) + present = {axis for dim, axis in dim_to_axis.items() if dim in ds_dims} + return sorted(present, key=_XGRID_AXES_ORDERING.index) +``` + +Replace `assert_all_field_dims_have_axis` (lines 54-76). Note: `model.py` imports this by name so the name must not change: + +```python +def assert_all_field_dims_have_axis(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata) -> None: + dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} + ax_dims = [(dim_to_axis.get(str(dim)), str(dim)) for dim in da.dims] + + for ax, dim_name in ax_dims: + if ax is None: + raise ValueError( + f'Dimension "{dim_name}" has no axis attribute. ' + f'HINT: You may want to add an {{"axis": A}} to your DataSet["{dim_name}"], where A is one of "X", "Y", "Z" or "T"' + ) + + seen_axes: dict[str, str] = {} + for ax, dim_name in ax_dims: + if ax in seen_axes: + raise ValueError( + f"Two dimensions ({dim_name!r} and {seen_axes[ax]!r}) provide values in the axis direction {ax!r}. " + "This is not possible, a field cannot have two dimensions on a single axis." + ) + seen_axes[ax] = dim_name + assert len(ax_dims) <= 4, ( + "The input dataset appears to have more than 4 dimensions after conversion. Execution should never reach this point. Please file an issue sharing more about your input dataset." + ) +``` + +- [ ] **Step 4: Delete the adapter classes and construct function** + +Delete `XgcmLikeAxis` (lines 116-118), `XgcmLikeGrid` (lines 121-129), and `construct_xgcm_axes_object` (lines 132-159) entirely. + +- [ ] **Step 5: Rewrite `XGrid.__init__`** + +Replace the `__init__` body from `grid = XgcmLikeGrid(...)` onward: + +```python +def __init__(self, model_data: xr.Dataset, mesh: Literal["flat", "spherical"] | SphericalMesh): + self.sgrid_metadata = model_data.sgrid.metadata + self._ds = model_data + self._mesh = get_mesh(mesh) + self._spatialhash = None + ds = model_data + + if "lon" in ds: + ds.set_coords("lon") + if "lat" in ds: + ds.set_coords("lat") + + axes = self.axes # uses _get_xgrid_axes(self.sgrid_metadata, set(self._ds.dims)) + if len(set(axes) & {"X", "Y"}) > 0: + assert_valid_lat_lon(ds["lat"], ds["lon"], self.sgrid_metadata) + + if "Z" in axes: + assert_valid_depth(ds["depth"]) + + self._ds = ds +``` + +- [ ] **Step 6: Rewrite `XGrid` properties and methods** + +Replace the `axes` property: + +```python +@property +def axes(self) -> list[ptyping.XgridAxis]: + return _get_xgrid_axes(self.sgrid_metadata, set(self._ds.dims)) +``` + +Replace the `lon` property: + +```python +@property +def lon(self): + """ + Note + ---- + Included for compatibility with v3 codebase. May be removed in future. + TODO v4: Evaluate + """ + if "X" not in self.axes: + return np.zeros(1) + if is_dask_collection(self._ds["lon"].data): + self._ds["lon"].load() + return self._ds["lon"].values +``` + +Replace the `lat` property: + +```python +@property +def lat(self): + """ + Note + ---- + Included for compatibility with v3 codebase. May be removed in future. + TODO v4: Evaluate + """ + if "Y" not in self.axes: + return np.zeros(1) + if is_dask_collection(self._ds["lat"].data): + self._ds["lat"].load() + return self._ds["lat"].values +``` + +Replace the `depth` property: + +```python +@property +def depth(self): + """ + Note + ---- + Included for compatibility with v3 codebase. May be removed in future. + TODO v4: Evaluate + """ + if "Z" not in self.axes: + return np.zeros(1) + return self._ds["depth"].values +``` + +Replace the `_datetimes` property: + +```python +@property +def _datetimes(self): + if "time" not in self._ds.dims: + return np.zeros(1) + return get_time(self._ds, "time") +``` + +Replace `get_axis_dim`: + +```python +def get_axis_dim(self, axis: ptyping.XgridAxis) -> int: + if axis not in self.axes: + raise ValueError(f"Axis {axis!r} is not part of this grid. Available axes: {self.axes}") + + fnp_x, fnp_y = self.sgrid_metadata.face_dimensions + if axis == "X": + return get_cell_count_along_dim(self._ds, fnp_x) + if axis == "Y": + return get_cell_count_along_dim(self._ds, fnp_y) + # axis == "Z" + assert self.sgrid_metadata.vertical_dimensions is not None + return get_cell_count_along_dim(self._ds, self.sgrid_metadata.vertical_dimensions[0]) +``` + +Replace `localize`: + +```python +def localize( + self, position: dict[ptyping.XgridAxis, tuple[int, float]], dims: list[str] +) -> dict[str, tuple[int, float]]: + """ + Uses the grid context (i.e., the staggering of the grid) to convert a position relative + to the F-points in the grid to a position relative to the staggered grid the array + of interest is defined on. + + Uses dimensions of the DataArray to determine the staggered grid. + + WARNING: This API is unstable and subject to change in future versions. + + Parameters + ---------- + position : dict + A mapping of the axis to a tuple of (index, barycentric coordinate) for the + F-points in the grid. + dims : list[str] + A list of dimension names that the DataArray is defined on. This is used to determine + the staggering of the grid and which axis each dimension corresponds to. + + Returns + ------- + dict[str, tuple[int, float]] + A mapping of the dimension names to a tuple of (index, barycentric coordinate) for + the staggered grid the DataArray is defined on. + + Example + ------- + >>> position = {'X': (5, 0.51), 'Y': ( + 10, 0.25), 'Z': (3, 0.75)} + >>> dims = ['time', 'depth', 'YC', 'XC'] + >>> grid.localize(position, dims) + {'depth': (3, 0.75), 'YC': (9, 0.75), 'XC': (5, 0.01)} + """ + dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) | {"time": "T"} + axis_to_var = {dim_to_axis[dim]: dim for dim in dims if dim in dim_to_axis} + var_positions = { + axis: get_dim_position(self.sgrid_metadata, dim) + for axis, dim in axis_to_var.items() + if axis != "T" + } + return { + axis_to_var[axis]: _convert_center_pos_to_fpoint( + index=index, + bcoord=bcoord, + position=var_positions[axis], + f_point_position=self._fpoint_info[axis], + ) + for axis, (index, bcoord) in position.items() + } +``` + +Replace `_fpoint_info`: + +```python +@cached_property +def _fpoint_info(self) -> dict[ptyping.XgridAxis, sgrid.Padding]: + """Returns a mapping of the spatial axes in the Grid to their Padding values (node positions).""" + metadata = self.sgrid_metadata + fnp_x, fnp_y = metadata.face_dimensions + result: dict[ptyping.XgridAxis, sgrid.Padding] = {} + axes = self.axes + if "X" in axes: + result["X"] = fnp_x.padding + if "Y" in axes: + result["Y"] = fnp_y.padding + if "Z" in axes and metadata.vertical_dimensions: + result["Z"] = metadata.vertical_dimensions[0].padding + return result +``` + +Replace `get_axis_dim_mapping`: + +```python +def get_axis_dim_mapping(self, dims: Sequence[Hashable]) -> dict[ptyping.XgridAxis, str]: + """ + Maps xarray dimension names to their corresponding axis (X, Y, Z). + + WARNING: This API is unstable and subject to change in future versions. + + Parameters + ---------- + dims : Sequence[Hashable] + Sequence of xarray dimension names + + Returns + ------- + dict[_XGRID_AXES, str] + Dictionary mapping axes (X, Y, Z) to their corresponding dimension names + + Examples + -------- + >>> grid.get_axis_dim_mapping(['time', 'lat', 'lon']) + {'Y': 'lat', 'X': 'lon'} + + Notes + ----- + Only returns mappings for spatial axes (X, Y, Z) that are present in the grid. + """ + dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) + result = {} + for dim in dims: + axis = dim_to_axis.get(str(dim)) + if axis in self.axes: + result[cast(ptyping.XgridAxis, axis)] = str(dim) + return result +``` + +- [ ] **Step 7: Rewrite module-level functions below `XGrid`** + +Replace `get_axis_from_dim_name` (lines 459-464): + +```python +def get_axis_from_dim_name(metadata: sgrid.SGrid2DMetadata, dim: Hashable) -> ptyping.XgcmAxisDirection | None: + """For a given dimension name in a grid, returns the direction axis it is on.""" + dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} + return dim_to_axis.get(str(dim)) +``` + +Replace `get_xgcm_position_from_dim_name` (lines 467-473). Keep the name as-is since it may be imported externally: + +```python +def get_xgcm_position_from_dim_name(metadata: sgrid.SGrid2DMetadata, dim: str) -> ptyping.GridPosition | None: + """For a given dimension, returns the GridPosition of the variable in the grid.""" + try: + return get_dim_position(metadata, dim) + except ValueError: + return None +``` + +Replace `assert_all_dimensions_correspond_with_axis` (lines 477-484): + +```python +def assert_all_dimensions_correspond_with_axis(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata) -> None: + dim_to_axis = _get_dim_to_axis_mapping(metadata) + for dim in da.dims: + if dim not in dim_to_axis: + raise ValueError( + f"Dimension {dim!r} for DataArray {da.name!r} with dims {da.dims} is not associated with a direction on the provided grid." + ) +``` + +Replace `assert_valid_field_array` (lines 487-509): + +```python +def assert_valid_field_array(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata): + """ + Asserts that for a data array: + - All dimensions are associated with a direction on the grid + - These directions are T, Z, Y, X and the array is ordered as T, Z, Y, X + """ + dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} + + for dim in da.dims: + if dim not in dim_to_axis: + raise ValueError( + f"Dimension {dim!r} for DataArray {da.name!r} with dims {da.dims} is not associated with a direction on the provided grid." + ) + + dim_to_axis_for_da = {dim: dim_to_axis[dim] for dim in da.dims} + dim_to_axis_for_da = cast(dict[Hashable, ptyping.XgcmAxisDirection], dim_to_axis_for_da) + + if set(dim_to_axis_for_da.values()) != {"T", "Z", "Y", "X"}: + raise ValueError( + f"DataArray {da.name!r} with dims {da.dims} has directions {tuple(dim_to_axis_for_da.values())}." + "Expected directions of 'T', 'Z', 'Y', and 'X'." + ) + + if list(dim_to_axis_for_da.values()) != ["T", "Z", "Y", "X"]: + raise ValueError( + f"Dimension order for array {da.name!r} is not valid. Got {tuple(dim_to_axis_for_da.keys())} with associated directions of {tuple(dim_to_axis_for_da.values())}. Expected directions of ('T', 'Z', 'Y', 'X'). Transpose your array accordingly." + ) +``` + +Replace `assert_valid_lat_lon` (lines 512-579): + +```python +def assert_valid_lat_lon(da_lat, da_lon, metadata: sgrid.SGrid2DMetadata): + """ + Asserts that the provided longitude and latitude DataArrays are defined appropriately + on the F points to match the internal representation in Parcels. + + - Longitude and latitude must be 1D or 2D (both must have the same dimensionality) + - Both are defined on the node points (i.e., not the face/center) + - If 1D: + - Longitude is associated with the X axis + - Latitude is associated with the Y axis + - If 2D: + - Lon and lat are defined on the same dimensions + - Lon and lat are transposed such they're Y, X + """ + assert_all_dimensions_correspond_with_axis(da_lon, metadata) + assert_all_dimensions_correspond_with_axis(da_lat, metadata) + + for dim in da_lon.dims: + if get_dim_position(metadata, dim) == "face": + raise ValueError( + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is defined on the center of the grid, but must be defined on the F points." + ) + for dim in da_lat.dims: + if get_dim_position(metadata, dim) == "face": + raise ValueError( + f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is defined on the center of the grid, but must be defined on the F points." + ) + + if da_lon.ndim != da_lat.ndim: + raise ValueError( + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} have different dimensionalities." + ) + if da_lon.ndim not in (1, 2): + raise ValueError( + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be 1D or 2D." + ) + + dim_to_axis = _get_dim_to_axis_mapping(metadata) + + if da_lon.ndim == 1: + if dim_to_axis.get(da_lon.dims[0]) != "X": + raise ValueError( + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is not associated with the X axis." + ) + if dim_to_axis.get(da_lat.dims[0]) != "Y": + raise ValueError( + f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is not associated with the Y axis." + ) + + if not np.all(np.diff(da_lon.values) > 0): + raise ValueError( + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} must be strictly increasing." + ) + if not np.all(np.diff(da_lat.values) > 0): + raise ValueError(f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be strictly increasing.") + + if da_lon.ndim == 2: + if da_lon.dims != da_lat.dims: + raise ValueError( + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be defined on the same dimensions." + ) + + lon_axes = [dim_to_axis.get(dim) for dim in da_lon.dims] + if lon_axes != ["Y", "X"]: + raise ValueError( + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be defined on the X and Y axes and transposed to have dimensions in order of Y, X." + ) +``` + +Replace `_convert_center_pos_to_fpoint` (lines 590-616): + +```python +def _convert_center_pos_to_fpoint( + *, + index: int, + bcoord: float, + position: ptyping.GridPosition, + f_point_position: sgrid.Padding, +) -> tuple[int, float]: + """Converts a physical position relative to the cell edges defined in the grid to be relative to the center point. + + This is used to "localize" a position to be relative to the staggered grid at which the field is defined, so that + it can be easily interpolated. + + This also handles different model input cell edges and centers are staggered in different directions (e.g., with NEMO and MITgcm). + """ + if position != "face": # Data is already defined on the F points + return index, bcoord + + bcoord = bcoord - 0.5 + if bcoord < 0: + bcoord += 1.0 + index -= 1 + + # Correct relative to the f-point position + # Padding.BOTH was "inner", Padding.LOW was "right" in xgcm vocabulary + if f_point_position in (sgrid.Padding.BOTH, sgrid.Padding.LOW): + index += 1 + + return index, bcoord +``` + +- [ ] **Step 8: Clean up `_typing.py`** + +Remove these lines from `src/parcels/_typing.py`: + +```python +# DELETE the TYPE_CHECKING block for xgcm: +if TYPE_CHECKING: + import xgcm + +# DELETE these two type aliases: +XgcmAxisPosition = Literal["center", "left", "right", "inner", "outer"] +XgcmAxes = Mapping[XgcmAxisDirection, "xgcm.Axis"] +``` + +Also remove `Mapping` from the `collections.abc` import if it is now unused (check by searching for other uses of `Mapping` in the file). + +- [ ] **Step 9: Run tests** + +```bash +pytest tests/sgrid/ tests/datasets/ -x -q +``` + +Expected: all pass. If any fail, fix before proceeding. + +- [ ] **Step 10: Commit** + +```bash +git add src/parcels/_core/xgrid.py src/parcels/_typing.py +git commit -m "refactor: replace XgcmLike adapter layer with direct SGRID metadata usage + +Co-authored-by: Claude " +``` + +--- + +### Task 3: Remove bridge functions, fix tests, remove xgcm dep + +Delete the now-dead bridge code in `_sgrid/core.py`, fix the two test files, and drop xgcm from `pyproject.toml`. + +**Files:** + +- Modify: `src/parcels/_sgrid/core.py` +- Modify: `src/parcels/_sgrid/__init__.py` +- Modify: `tests/sgrid/test_sgrid.py` +- Modify: `tests/datasets/test_structured.py` +- Modify: `pyproject.toml` + +**Interfaces:** + +- Consumes: `get_dim_position` from Task 1 + +- [ ] **Step 1: Remove `SGRID_PADDING_TO_XGCM_POSITION` and `xgcm_parse_sgrid` from `_sgrid/core.py`** + +Delete lines 41-47 (the `SGRID_PADDING_TO_XGCM_POSITION` dict): + +```python +# DELETE: +SGRID_PADDING_TO_XGCM_POSITION = { + Padding.LOW: "right", + Padding.HIGH: "left", + Padding.BOTH: "inner", + Padding.NONE: "outer", + # "center" position is not used in SGrid, in SGrid this would just be the edges/faces themselves +} +``` + +Delete lines 470-492 (the `xgcm_parse_sgrid` function): + +```python +# DELETE: +def xgcm_parse_sgrid(ds: xr.Dataset): + # Function similar to that provided in `xgcm.metadata_parsers. + # Might at some point be upstreamed to xgcm directly + grid = ds.sgrid.metadata + ... + return (ds, {"coords": xgcm_coords}) +``` + +- [ ] **Step 2: Update `_sgrid/__init__.py`** + +Remove `xgcm_parse_sgrid` from the import and `__all__`: + +```python +# Change from: +from .core import ( + FaceNodePadding, + Padding, + SGrid2DMetadata, + SGrid3DMetadata, + _attach_sgrid_metadata, + dump_mappings, + get_n_faces, + get_n_nodes, + load_mappings, + xgcm_parse_sgrid, +) + +__all__ = [ + "FaceNodePadding", + "Padding", + "SGrid2DMetadata", + "SGrid3DMetadata", + "SgridAccessor", + "_attach_sgrid_metadata", + "dump_mappings", + "get_n_faces", + "get_n_nodes", + "load_mappings", + "xgcm_parse_sgrid", +] + +# To: +from .core import ( + FaceNodePadding, + Padding, + SGrid2DMetadata, + SGrid3DMetadata, + _attach_sgrid_metadata, + dump_mappings, + get_n_faces, + get_n_nodes, + load_mappings, +) + +__all__ = [ + "FaceNodePadding", + "Padding", + "SGrid2DMetadata", + "SGrid3DMetadata", + "SgridAccessor", + "_attach_sgrid_metadata", + "dump_mappings", + "get_n_faces", + "get_n_nodes", + "load_mappings", +] +``` + +- [ ] **Step 3: Fix `tests/sgrid/test_sgrid.py`** + +Remove the `import xgcm` line (line 7). + +Remove `SGRID_PADDING_TO_XGCM_POSITION` from the import on line 12: + +```python +# Change from: +from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION, _get_unique_names, parse_grid_attrs + +# To: +from parcels._sgrid.core import _get_unique_names, parse_grid_attrs +``` + +Delete the entire `test_parse_sgrid_2d` function (lines 253-273) and `test_parse_sgrid_3d` function (lines 276-290) — these test the now-deleted `xgcm_parse_sgrid` bridge function. The SGRID parsing itself is still exercised by the existing `test_Grid2DMetadata_roundtrip`, `test_parse_grid_attrs`, and other tests. + +- [ ] **Step 4: Fix `tests/datasets/test_structured.py`** + +Replace the entire file content: + +```python +import parcels._sgrid as sgrid +from parcels._datasets.structured.generic import datasets +from parcels._sgrid.accessor import get_dim_position + + +def test_left_indexed_dataset(): + """Checks that 'ds_2d_left' uses HIGH padding (MITgcm / left-indexed style).""" + ds = datasets["ds_2d_left"] + metadata = ds.sgrid.metadata + for fnp in metadata.face_dimensions: + assert get_dim_position(metadata, fnp.face) == "face" + assert get_dim_position(metadata, fnp.node) == sgrid.Padding.HIGH + + +def test_right_indexed_dataset(): + """Checks that 'ds_2d_right' uses LOW padding (NEMO / right-indexed style).""" + ds = datasets["ds_2d_right"] + metadata = ds.sgrid.metadata + for fnp in metadata.face_dimensions: + assert get_dim_position(metadata, fnp.face) == "face" + assert get_dim_position(metadata, fnp.node) == sgrid.Padding.LOW +``` + +- [ ] **Step 5: Remove xgcm from `pyproject.toml`** + +Delete the line `"xgcm >=0.9.0",` from the `dependencies` list in `pyproject.toml`. + +- [ ] **Step 6: Run the full test subset** + +```bash +pytest tests/sgrid/ tests/datasets/ -x -q +``` + +Expected: all pass. The two deleted tests (`test_parse_sgrid_2d`, `test_parse_sgrid_3d`) are gone; four new tests from Task 1 and two rewritten tests from `test_structured.py` all pass. + +- [ ] **Step 7: Verify xgcm is no longer imported anywhere in src/** + +```bash +grep -r "import xgcm" src/parcels --include="*.py" +``` + +Expected: no output. + +- [ ] **Step 8: Commit** + +```bash +git add src/parcels/_sgrid/core.py src/parcels/_sgrid/__init__.py \ + tests/sgrid/test_sgrid.py tests/datasets/test_structured.py \ + pyproject.toml +git commit -m "feat: remove xgcm dependency — use SGRID metadata natively throughout + +Co-authored-by: Claude " +``` From 07550f1431bbdbb81498424c051b5d435a3f8a5c Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:33:38 +0800 Subject: [PATCH 03/13] feat: add GridPosition type and get_dim_position() SGRID helper --- src/parcels/_sgrid/accessor.py | 13 +++++++++++++ src/parcels/_typing.py | 2 ++ tests/sgrid/test_sgrid.py | 31 +++++++++++++++++++++++++++++++ 3 files changed, 46 insertions(+) diff --git a/src/parcels/_sgrid/accessor.py b/src/parcels/_sgrid/accessor.py index a65f52e4d..62fede102 100644 --- a/src/parcels/_sgrid/accessor.py +++ b/src/parcels/_sgrid/accessor.py @@ -148,6 +148,19 @@ def _get_axis_info(grid: SGrid2DMetadata) -> dict[str, tuple[FaceNodePadding, bo return result +def get_dim_position(grid: SGrid2DMetadata, dim: str) -> "Literal['face'] | Padding": + """Returns 'face' if dim is a face dimension, or the Padding value if it is a node dimension. + + Replaces xgcm's position string vocabulary ('center', 'left', 'right', 'inner', 'outer') + with SGRID-native types. + """ + axis_info = _get_axis_info(grid) + if dim not in axis_info: + raise ValueError(f"Dimension {dim!r} is not a spatial SGRID dimension in this grid.") + fnp, is_node = axis_info[dim] + return fnp.padding if is_node else "face" + + def _derive_paired_indexer( indexer: Any, indexer_is_node: bool, diff --git a/src/parcels/_typing.py b/src/parcels/_typing.py index 87e23a53c..24c10502c 100644 --- a/src/parcels/_typing.py +++ b/src/parcels/_typing.py @@ -15,6 +15,7 @@ from cftime import datetime as cftime_datetime from parcels._core.mesh import TMesh # noqa: F401 +from parcels._sgrid.core import Padding if TYPE_CHECKING: import xgcm @@ -48,6 +49,7 @@ CfAxis = XgcmAxisDirection XgcmAxisPosition = Literal["center", "left", "right", "inner", "outer"] XgcmAxes = Mapping[XgcmAxisDirection, "xgcm.Axis"] +GridPosition = Literal["face"] | Padding VectorFields = dict[str, tuple[str, str] | tuple[str, str, str]] diff --git a/tests/sgrid/test_sgrid.py b/tests/sgrid/test_sgrid.py index 5582821d0..b83bff494 100644 --- a/tests/sgrid/test_sgrid.py +++ b/tests/sgrid/test_sgrid.py @@ -9,6 +9,7 @@ import parcels._sgrid as sgrid import parcels._strategies as pst +from parcels._sgrid.accessor import get_dim_position from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION, _get_unique_names, parse_grid_attrs @@ -78,6 +79,36 @@ def test_get_value_by_id(sgrid_metadata: sgrid.SGrid2DMetadata | sgrid.SGrid3DMe assert sgrid_metadata.get_value_by_id(id) == value +def test_get_dim_position_face_dims(): + """Face dimensions return 'face'.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) + # face_dimensions = (FaceNodePadding("face_dimension1", "node_dimension1", Padding.LOW), ...) + assert get_dim_position(metadata, "face_dimension1") == "face" + assert get_dim_position(metadata, "face_dimension2") == "face" + + +def test_get_dim_position_node_dims(): + """Node dimensions return their Padding value.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) + assert get_dim_position(metadata, "node_dimension1") == sgrid.Padding.LOW + assert get_dim_position(metadata, "node_dimension2") == sgrid.Padding.LOW + + +def test_get_dim_position_vertical(): + """Vertical face and node dimensions are handled.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=True, with_node_coordinates=False) + # vertical_dimensions = (FaceNodePadding("vertical_dimensions_dim1", "vertical_dimensions_dim2", Padding.LOW),) + assert get_dim_position(metadata, "vertical_dimensions_dim1") == "face" + assert get_dim_position(metadata, "vertical_dimensions_dim2") == sgrid.Padding.LOW + + +def test_get_dim_position_unknown_dim_raises(): + """Unknown dimensions raise ValueError.""" + metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) + with pytest.raises(ValueError, match="not a spatial SGRID dimension"): + get_dim_position(metadata, "nonexistent_dim") + + def dummy_sgrid_ds(grid: sgrid.SGrid2DMetadata | sgrid.SGrid3DMetadata) -> xr.Dataset: if isinstance(grid, sgrid.SGrid2DMetadata): return dummy_sgrid_2d_ds(grid) From d6b8d8ccaa0b57956e94696f79d24e3ed09df1a5 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:39:02 +0800 Subject: [PATCH 04/13] refactor: replace XgcmLike adapter layer with direct SGRID metadata usage Remove XgcmLikeAxis, XgcmLikeGrid, construct_xgcm_axes_object and all xgcm imports from xgrid.py; replace with SGrid2DMetadata / FaceNodePadding / Padding throughout. Clean up XgcmAxisPosition and XgcmAxes from _typing.py. Update test_structured.py to use SGRID-native get_dim_position assertions. --- src/parcels/_core/xgrid.py | 259 ++++++++++++------------------ src/parcels/_typing.py | 9 +- tests/datasets/test_structured.py | 26 ++- 3 files changed, 117 insertions(+), 177 deletions(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 5f76ea9e4..345c2a211 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -1,12 +1,10 @@ from collections.abc import Hashable, Sequence -from dataclasses import dataclass from functools import cached_property -from typing import TYPE_CHECKING, Any, Literal, cast +from typing import Literal, cast import numpy as np import numpy.typing as npt import xarray as xr -import xgcm from dask import is_dask_collection import parcels._sgrid as sgrid @@ -14,32 +12,24 @@ from parcels._core.basegrid import BaseGrid from parcels._core.index_search import _search_1d_array, _search_indices_curvilinear_2d from parcels._core.mesh import SphericalMesh, get_mesh -from parcels._sgrid.accessor import _get_dim_to_axis_mapping -from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION - -if TYPE_CHECKING: - import xgcm.axis +from parcels._sgrid.accessor import _get_dim_to_axis_mapping, get_dim_position _FIELD_DATA_ORDERING: Sequence[ptyping.XgcmAxisDirection] = "TZYX" _XGRID_AXES_ORDERING: Sequence[ptyping.XgridAxis] = "ZYX" -_DEFAULT_XGCM_KWARGS: dict[str, Any] = {"padding": "fill"} - -def get_cell_count_along_dim(ds: xr.Dataset, axis: xgcm.axis.Axis) -> int: - first_coord = list(axis.coords.items())[0] - _, coord_var = first_coord +def get_cell_count_along_dim(ds: xr.Dataset, fnp: sgrid.FaceNodePadding) -> int: + return ds[fnp.face].size - 1 - return ds[coord_var].size - 1 +def get_time(ds: xr.Dataset, time_dim: str) -> npt.NDArray: + return ds[time_dim].values -def get_time(ds: xr.Dataset, axis: xgcm.axis.Axis) -> npt.NDArray: - return ds[axis.coords["center"]].values - -def _get_xgrid_axes(grid: xgcm.Grid) -> list[ptyping.XgridAxis]: - spatial_axes = [a for a in grid.axes.keys() if a in ["X", "Y", "Z"]] - return sorted(spatial_axes, key=_XGRID_AXES_ORDERING.index) +def _get_xgrid_axes(metadata: sgrid.SGrid2DMetadata, ds_dims: set[str]) -> list[ptyping.XgridAxis]: + dim_to_axis = _get_dim_to_axis_mapping(metadata) + present = {axis for dim, axis in dim_to_axis.items() if dim in ds_dims} + return sorted(present, key=_XGRID_AXES_ORDERING.index) def _drop_field_data(ds: xr.Dataset) -> xr.Dataset: @@ -51,17 +41,17 @@ def _drop_field_data(ds: xr.Dataset) -> xr.Dataset: return ds.drop_vars(set(ds.data_vars) - {"grid"}) # don't drop sgrid metadata -def assert_all_field_dims_have_axis(da: xr.DataArray, xgcm_grid: xgcm.Grid) -> None: - ax_dims = [(get_axis_from_dim_name(xgcm_grid.axes, dim), dim) for dim in da.dims] - for dim in ax_dims: - if dim[0] is None: +def assert_all_field_dims_have_axis(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata) -> None: + dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} + ax_dims = [(dim_to_axis.get(str(dim)), str(dim)) for dim in da.dims] + + for ax, dim_name in ax_dims: + if ax is None: raise ValueError( - f'Dimension "{dim[1]}" has no axis attribute. ' - f'HINT: You may want to add an {{"axis": A}} to your DataSet["{dim[1]}"], where A is one of "X", "Y", "Z" or "T"' + f'Dimension "{dim_name}" has no axis attribute. ' + f'HINT: You may want to add an {{"axis": A}} to your DataSet["{dim_name}"], where A is one of "X", "Y", "Z" or "T"' ) - ax_dims = cast(dict[str, str], ax_dims) - seen_axes: dict[str, str] = {} for ax, dim_name in ax_dims: if ax in seen_axes: @@ -113,57 +103,12 @@ def _transpose_xfield_data_to_tzyx(da: xr.DataArray, sgrid_metadata: sgrid.SGrid return da.transpose(*[ax_dim[0] for ax_dim in ax_dims]) -@dataclass -class XgcmLikeAxis: - coords: dict[ptyping.XgcmAxisPosition, str] - - -class XgcmLikeGrid: - """Adapter class to circumvent XGCM as a dep. - - - TODO: This is only used as a temporary class for the moment. Down the line we should refactor to remove XGCM entirely and work with SGRID metadata (especially since COMODO metadata isn't standard, and nor is the xgcm data model). - """ - - def __init__(self, sgrid_metadata: sgrid.SGrid2DMetadata, model_data: xr.Dataset): - self.axes: dict[ptyping.CfAxisSpatial, XgcmLikeAxis] = construct_xgcm_axes_object(sgrid_metadata, model_data) - - -def construct_xgcm_axes_object(metadata: sgrid.SGrid2DMetadata, model_data: xr.Dataset) -> dict[str, XgcmLikeAxis]: - lst: list[tuple[ptyping.CfAxis, str, ptyping.XgcmAxisPosition]] = [] - - for fnp, axis in zip(metadata.face_dimensions, ("X", "Y"), strict=True): - lst.append((axis, fnp.face, "center")) - lst.append((axis, fnp.node, SGRID_PADDING_TO_XGCM_POSITION[fnp.padding])) - - if metadata.vertical_dimensions is not None: - assert len(metadata.vertical_dimensions) == 1 - fnp = metadata.vertical_dimensions[0] - axis = "Z" - lst.append((axis, fnp.face, "center")) - lst.append((axis, fnp.node, SGRID_PADDING_TO_XGCM_POSITION[fnp.padding])) - - # filter so that only dims in the dataset itself are mentioned - lst = [i for i in lst if i[1] in model_data.dims] - - # Add time axis to xgcm_kwargs if present - if "time" in model_data.dims: - lst.append(("T", "time", "center")) - - ret = {} - for axis, dim, position in lst: - if axis not in ret: - ret[axis] = XgcmLikeAxis({}) - ret[axis].coords[position] = dim - - return ret - - class XGrid(BaseGrid): """ - Class to represent a structured grid in Parcels. Wraps a xgcm-like Grid object (we use a trimmed down version of the xgcm.Grid class that is vendored with Parcels). + Class to represent a structured grid in Parcels. This class provides methods and properties required for indexing and interpolating on the grid. + Grid topology is derived directly from SGRID metadata attached to the dataset. Assumptions: - If using Parcels in the context of a spatially periodic simulation, the provided grid already has a halo @@ -173,8 +118,6 @@ class XGrid(BaseGrid): def __init__(self, model_data: xr.Dataset, mesh: Literal["flat", "spherical"] | SphericalMesh): self.sgrid_metadata = model_data.sgrid.metadata self._ds = model_data - grid = XgcmLikeGrid(self.sgrid_metadata, model_data) - self.xgcm_grid = grid self._mesh = get_mesh(mesh) self._spatialhash = None ds = model_data @@ -185,10 +128,11 @@ def __init__(self, model_data: xr.Dataset, mesh: Literal["flat", "spherical"] | if "lat" in ds: ds.set_coords("lat") - if len(set(grid.axes) & {"X", "Y"}) > 0: # Only if spatial grid is >0D (see #2054 for further development) - assert_valid_lat_lon(ds["lat"], ds["lon"], grid.axes) + axes = self.axes + if len(set(axes) & {"X", "Y"}) > 0: # Only if spatial grid is >0D (see #2054 for further development) + assert_valid_lat_lon(ds["lat"], ds["lon"], self.sgrid_metadata) - if "Z" in grid.axes: + if "Z" in axes: assert_valid_depth(ds["depth"]) self._ds = ds @@ -198,7 +142,7 @@ def __init__(self, model_data: xr.Dataset, mesh: Literal["flat", "spherical"] | @property def axes(self) -> list[ptyping.XgridAxis]: - return _get_xgrid_axes(self.xgcm_grid) + return _get_xgrid_axes(self.sgrid_metadata, set(self._ds.dims)) @property def lon(self): @@ -208,9 +152,7 @@ def lon(self): Included for compatibility with v3 codebase. May be removed in future. TODO v4: Evaluate """ - try: - _ = self.xgcm_grid.axes["X"] - except KeyError: + if "X" not in self.axes: return np.zeros(1) # ensure lon is loaded into memory for dask-backed datasets, as it is used in the search method if is_dask_collection(self._ds["lon"].data): @@ -225,9 +167,7 @@ def lat(self): Included for compatibility with v3 codebase. May be removed in future. TODO v4: Evaluate """ - try: - _ = self.xgcm_grid.axes["Y"] - except KeyError: + if "Y" not in self.axes: return np.zeros(1) # ensure lat is loaded into memory for dask-backed datasets, as it is used in the search method if is_dask_collection(self._ds["lat"].data): @@ -242,19 +182,15 @@ def depth(self): Included for compatibility with v3 codebase. May be removed in future. TODO v4: Evaluate """ - try: - _ = self.xgcm_grid.axes["Z"] - except KeyError: + if "Z" not in self.axes: return np.zeros(1) return self._ds["depth"].values @property def _datetimes(self): - try: - axis = self.xgcm_grid.axes["T"] - except KeyError: + if "time" not in self._ds.dims: return np.zeros(1) - return get_time(self._ds, axis) + return get_time(self._ds, "time") @property def time(self): @@ -283,7 +219,14 @@ def get_axis_dim(self, axis: ptyping.XgridAxis) -> int: if axis not in self.axes: raise ValueError(f"Axis {axis!r} is not part of this grid. Available axes: {self.axes}") - return get_cell_count_along_dim(self._ds, self.xgcm_grid.axes[axis]) + fnp_x, fnp_y = self.sgrid_metadata.face_dimensions + if axis == "X": + return get_cell_count_along_dim(self._ds, fnp_x) + if axis == "Y": + return get_cell_count_along_dim(self._ds, fnp_y) + # axis == "Z" + assert self.sgrid_metadata.vertical_dimensions is not None + return get_cell_count_along_dim(self._ds, self.sgrid_metadata.vertical_dimensions[0]) def localize( self, position: dict[ptyping.XgridAxis, tuple[int, float]], dims: list[str] @@ -320,16 +263,17 @@ def localize( >>> grid.localize(position, dims) {'depth': (3, 0.75), 'YC': (9, 0.75), 'XC': (5, 0.01)} """ - axis_to_var = {get_axis_from_dim_name(self.xgcm_grid.axes, dim): dim for dim in dims} + dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) | {"time": "T"} + axis_to_var = {dim_to_axis[dim]: dim for dim in dims if dim in dim_to_axis} var_positions = { - axis: get_xgcm_position_from_dim_name(self.xgcm_grid.axes, dim) for axis, dim in axis_to_var.items() + axis: get_dim_position(self.sgrid_metadata, dim) for axis, dim in axis_to_var.items() if axis != "T" } return { axis_to_var[axis]: _convert_center_pos_to_fpoint( index=index, bcoord=bcoord, - xgcm_position=var_positions[axis], - f_points_xgcm_position=self._fpoint_info[axis], + position=var_positions[axis], + f_point_position=self._fpoint_info[axis], ) for axis, (index, bcoord) in position.items() } @@ -410,18 +354,19 @@ def search(self, z, y, x, ei=None): } @cached_property - def _fpoint_info(self): - """Returns a mapping of the spatial axes in the Grid to their XGCM positions.""" - xgcm_axes = self.xgcm_grid.axes - f_point_positions = ["left", "right", "inner", "outer"] - axis_position_mapping = {} - for axis in self.axes: - coords = xgcm_axes[axis].coords - edge_positions = [pos for pos in coords.keys() if pos in f_point_positions] - assert len(edge_positions) == 1, f"Axis {axis} has multiple edge positions: {edge_positions}" - axis_position_mapping[axis] = edge_positions[0] - - return axis_position_mapping + def _fpoint_info(self) -> dict[ptyping.XgridAxis, sgrid.Padding]: + """Returns a mapping of the spatial axes in the Grid to their Padding values (node positions).""" + metadata = self.sgrid_metadata + fnp_x, fnp_y = metadata.face_dimensions + result: dict[ptyping.XgridAxis, sgrid.Padding] = {} + axes = self.axes + if "X" in axes: + result["X"] = fnp_x.padding + if "Y" in axes: + result["Y"] = fnp_y.padding + if "Z" in axes and metadata.vertical_dimensions: + result["Z"] = metadata.vertical_dimensions[0].padding + return result def get_axis_dim_mapping(self, dims: Sequence[Hashable]) -> dict[ptyping.XgridAxis, str]: """ @@ -448,74 +393,76 @@ def get_axis_dim_mapping(self, dims: Sequence[Hashable]) -> dict[ptyping.XgridAx ----- Only returns mappings for spatial axes (X, Y, Z) that are present in the grid. """ + dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) result = {} for dim in dims: - axis = get_axis_from_dim_name(self.xgcm_grid.axes, dim) - if axis in self.axes: # Only include spatial axes (X, Y, Z) - result[cast(ptyping.XgridAxis, axis)] = dim + axis = dim_to_axis.get(str(dim)) + if axis in self.axes: + result[cast(ptyping.XgridAxis, axis)] = str(dim) return result -def get_axis_from_dim_name(axes: ptyping.XgcmAxes, dim: Hashable) -> ptyping.XgcmAxisDirection | None: +def get_axis_from_dim_name(metadata: sgrid.SGrid2DMetadata, dim: Hashable) -> ptyping.XgcmAxisDirection | None: """For a given dimension name in a grid, returns the direction axis it is on.""" - for axis_name, axis in axes.items(): - if dim in axis.coords.values(): - return axis_name - return None - + dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} + return dim_to_axis.get(str(dim)) -def get_xgcm_position_from_dim_name(axes: ptyping.XgcmAxes, dim: str) -> ptyping.XgcmAxisPosition | None: - """For a given dimension, returns the position of the variable in the grid.""" - for axis in axes.values(): - var_to_position = {var: position for position, var in axis.coords.items()} - if dim in var_to_position: - return var_to_position[dim] # type: ignore[invalid-return-type] # due to mistyping in xgcm - return None +def get_xgcm_position_from_dim_name(metadata: sgrid.SGrid2DMetadata, dim: str) -> ptyping.GridPosition | None: + """For a given dimension, returns the GridPosition of the variable in the grid.""" + try: + return get_dim_position(metadata, dim) + except ValueError: + return None -def assert_all_dimensions_correspond_with_axis(da: xr.DataArray, axes: ptyping.XgcmAxes) -> None: - dim_to_axis = {dim: get_axis_from_dim_name(axes, dim) for dim in da.dims} - - for dim, direction in dim_to_axis.items(): - if direction is None: +def assert_all_dimensions_correspond_with_axis(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata) -> None: + dim_to_axis = _get_dim_to_axis_mapping(metadata) + for dim in da.dims: + if dim not in dim_to_axis: raise ValueError( f"Dimension {dim!r} for DataArray {da.name!r} with dims {da.dims} is not associated with a direction on the provided grid." ) -def assert_valid_field_array(da: xr.DataArray, axes: ptyping.XgcmAxes): +def assert_valid_field_array(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata): """ Asserts that for a data array: - All dimensions are associated with a direction on the grid - These directions are T, Z, Y, X and the array is ordered as T, Z, Y, X """ - assert_all_dimensions_correspond_with_axis(da, axes) + dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} - dim_to_axis = {dim: get_axis_from_dim_name(axes, dim) for dim in da.dims} - dim_to_axis = cast(dict[Hashable, ptyping.XgcmAxisDirection], dim_to_axis) + for dim in da.dims: + if dim not in dim_to_axis: + raise ValueError( + f"Dimension {dim!r} for DataArray {da.name!r} with dims {da.dims} is not associated with a direction on the provided grid." + ) + + dim_to_axis_for_da = {dim: dim_to_axis[dim] for dim in da.dims} + dim_to_axis_for_da = cast(dict[Hashable, ptyping.XgcmAxisDirection], dim_to_axis_for_da) # Assert all dimensions are present - if set(dim_to_axis.values()) != {"T", "Z", "Y", "X"}: + if set(dim_to_axis_for_da.values()) != {"T", "Z", "Y", "X"}: raise ValueError( - f"DataArray {da.name!r} with dims {da.dims} has directions {tuple(dim_to_axis.values())}." + f"DataArray {da.name!r} with dims {da.dims} has directions {tuple(dim_to_axis_for_da.values())}." "Expected directions of 'T', 'Z', 'Y', and 'X'." ) # Assert order is t, z, y, x - if list(dim_to_axis.values()) != ["T", "Z", "Y", "X"]: + if list(dim_to_axis_for_da.values()) != ["T", "Z", "Y", "X"]: raise ValueError( - f"Dimension order for array {da.name!r} is not valid. Got {tuple(dim_to_axis.keys())} with associated directions of {tuple(dim_to_axis.values())}. Expected directions of ('T', 'Z', 'Y', 'X'). Transpose your array accordingly." + f"Dimension order for array {da.name!r} is not valid. Got {tuple(dim_to_axis_for_da.keys())} with associated directions of {tuple(dim_to_axis_for_da.values())}. Expected directions of ('T', 'Z', 'Y', 'X'). Transpose your array accordingly." ) -def assert_valid_lat_lon(da_lat, da_lon, axes: ptyping.XgcmAxes): +def assert_valid_lat_lon(da_lat, da_lon, metadata: sgrid.SGrid2DMetadata): """ Asserts that the provided longitude and latitude DataArrays are defined appropriately on the F points to match the internal representation in Parcels. - Longitude and latitude must be 1D or 2D (both must have the same dimensionality) - - Both are defined on the left points (i.e., not the centers) + - Both are defined on the node points (i.e., not the face/center) - If 1D: - Longitude is associated with the X axis - Latitude is associated with the Y axis @@ -523,19 +470,16 @@ def assert_valid_lat_lon(da_lat, da_lon, axes: ptyping.XgcmAxes): - Lon and lat are defined on the same dimensions - Lon and lat are transposed such they're Y, X """ - assert_all_dimensions_correspond_with_axis(da_lon, axes) - assert_all_dimensions_correspond_with_axis(da_lat, axes) - - dim_to_position = {dim: get_xgcm_position_from_dim_name(axes, dim) for dim in da_lon.dims} - dim_to_position.update({dim: get_xgcm_position_from_dim_name(axes, dim) for dim in da_lat.dims}) + assert_all_dimensions_correspond_with_axis(da_lon, metadata) + assert_all_dimensions_correspond_with_axis(da_lat, metadata) for dim in da_lon.dims: - if get_xgcm_position_from_dim_name(axes, dim) == "center": + if get_dim_position(metadata, dim) == "face": raise ValueError( f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is defined on the center of the grid, but must be defined on the F points." ) for dim in da_lat.dims: - if get_xgcm_position_from_dim_name(axes, dim) == "center": + if get_dim_position(metadata, dim) == "face": raise ValueError( f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is defined on the center of the grid, but must be defined on the F points." ) @@ -549,12 +493,14 @@ def assert_valid_lat_lon(da_lat, da_lon, axes: ptyping.XgcmAxes): f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be 1D or 2D." ) + dim_to_axis = _get_dim_to_axis_mapping(metadata) + if da_lon.ndim == 1: - if get_axis_from_dim_name(axes, da_lon.dims[0]) != "X": + if dim_to_axis.get(da_lon.dims[0]) != "X": raise ValueError( f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is not associated with the X axis." ) - if get_axis_from_dim_name(axes, da_lat.dims[0]) != "Y": + if dim_to_axis.get(da_lat.dims[0]) != "Y": raise ValueError( f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is not associated with the Y axis." ) @@ -572,7 +518,7 @@ def assert_valid_lat_lon(da_lat, da_lon, axes: ptyping.XgcmAxes): f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be defined on the same dimensions." ) - lon_axes = [get_axis_from_dim_name(axes, dim) for dim in da_lon.dims] + lon_axes = [dim_to_axis.get(dim) for dim in da_lon.dims] if lon_axes != ["Y", "X"]: raise ValueError( f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be defined on the X and Y axes and transposed to have dimensions in order of Y, X." @@ -591,8 +537,8 @@ def _convert_center_pos_to_fpoint( *, index: int, bcoord: float, - xgcm_position: ptyping.XgcmAxisPosition, - f_points_xgcm_position: ptyping.XgcmAxisPosition, + position: ptyping.GridPosition, + f_point_position: sgrid.Padding, ) -> tuple[int, float]: """Converts a physical position relative to the cell edges defined in the grid to be relative to the center point. @@ -601,7 +547,7 @@ def _convert_center_pos_to_fpoint( This also handles different model input cell edges and centers are staggered in different directions (e.g., with NEMO and MITgcm). """ - if xgcm_position != "center": # Data is already defined on the F points + if position != "face": # Data is already defined on the F points return index, bcoord bcoord = bcoord - 0.5 @@ -610,7 +556,8 @@ def _convert_center_pos_to_fpoint( index -= 1 # Correct relative to the f-point position - if f_points_xgcm_position in ["inner", "right"]: + # Padding.BOTH was "inner", Padding.LOW was "right" in xgcm vocabulary + if f_point_position in (sgrid.Padding.BOTH, sgrid.Padding.LOW): index += 1 return index, bcoord diff --git a/src/parcels/_typing.py b/src/parcels/_typing.py index 24c10502c..dd98afc44 100644 --- a/src/parcels/_typing.py +++ b/src/parcels/_typing.py @@ -7,9 +7,9 @@ """ import os -from collections.abc import Callable, Mapping +from collections.abc import Callable from datetime import datetime -from typing import TYPE_CHECKING, Literal, get_args +from typing import Literal, get_args import numpy as np from cftime import datetime as cftime_datetime @@ -17,9 +17,6 @@ from parcels._core.mesh import TMesh # noqa: F401 from parcels._sgrid.core import Padding -if TYPE_CHECKING: - import xgcm - InterpMethodOption = Literal[ "linear", "nearest", @@ -47,8 +44,6 @@ XgridAxis = CfAxisSpatial XgcmAxisDirection = CfAxisSpatial | Literal["T"] CfAxis = XgcmAxisDirection -XgcmAxisPosition = Literal["center", "left", "right", "inner", "outer"] -XgcmAxes = Mapping[XgcmAxisDirection, "xgcm.Axis"] GridPosition = Literal["face"] | Padding VectorFields = dict[str, tuple[str, str] | tuple[str, str, str]] diff --git a/tests/datasets/test_structured.py b/tests/datasets/test_structured.py index 2750583f5..f96bd1636 100644 --- a/tests/datasets/test_structured.py +++ b/tests/datasets/test_structured.py @@ -1,23 +1,21 @@ -import xgcm - -from parcels._core.xgrid import _DEFAULT_XGCM_KWARGS from parcels._datasets.structured.generic import datasets +from parcels._sgrid.accessor import get_dim_position +from parcels._sgrid.core import Padding def test_left_indexed_dataset(): - """Checks that 'ds_2d_left' is right indexed on all variables.""" + """Checks that 'ds_2d_left' has HIGH padding (MITgcm/left-indexed) on all spatial axes.""" ds = datasets["ds_2d_left"] - grid = xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS) - - for _axis_name, axis in grid.axes.items(): - for pos, _dim_name in axis.coords.items(): - assert pos in ["left", "center"] + metadata = ds.sgrid.metadata + for fnp in metadata.face_dimensions: + assert get_dim_position(metadata, fnp.face) == "face" + assert get_dim_position(metadata, fnp.node) == Padding.HIGH def test_right_indexed_dataset(): - """Checks that 'ds_2d_right' is right indexed on all variables.""" + """Checks that 'ds_2d_right' has LOW padding (NEMO/right-indexed) on all spatial axes.""" ds = datasets["ds_2d_right"] - grid = xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS) - for _axis_name, axis in grid.axes.items(): - for pos, _dim_name in axis.coords.items(): - assert pos in ["center", "right"] + metadata = ds.sgrid.metadata + for fnp in metadata.face_dimensions: + assert get_dim_position(metadata, fnp.face) == "face" + assert get_dim_position(metadata, fnp.node) == Padding.LOW From 30113f13b559fe229927ec69738760cfeb6e4950 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:41:34 +0800 Subject: [PATCH 05/13] refactor: remove xgcm bridge functions and dependency Delete SGRID_PADDING_TO_XGCM_POSITION and xgcm_parse_sgrid() from _sgrid/core.py, remove xgcm_parse_sgrid from the package __all__, delete the two xgcm-backed tests from test_sgrid.py, and remove "xgcm >=0.9.0" from pyproject.toml and pixi.toml. xgcm is no longer referenced anywhere in src/parcels. --- pixi.toml | 2 -- pyproject.toml | 2 -- src/parcels/_sgrid/__init__.py | 2 -- src/parcels/_sgrid/core.py | 34 --------------------------- tests/sgrid/test_sgrid.py | 43 +--------------------------------- 5 files changed, 1 insertion(+), 82 deletions(-) diff --git a/pixi.toml b/pixi.toml index 50ca614ce..a3149893d 100644 --- a/pixi.toml +++ b/pixi.toml @@ -33,7 +33,6 @@ holoviews = ">=1.22.0" # https://github.com/prefix-dev/rattler-build/issues/2326 uxarray = ">=2026.04.1" dask = ">=2024.5.1" zarr = ">=3" -xgcm = { git = "https://github.com/VeckoTheGecko/xgcm", rev = "relax-pixi-build-python" } # TODO: Switch to release version after release is cut cf_xarray = ">=0.8.6" cftime = ">=1.6.3" pooch = ">=1.8.0" @@ -64,7 +63,6 @@ pyarrow = "20.0.*" uxarray = "==2026.04.1" dask = "2024.6.*" zarr = "3.0.*" -xgcm = { git = "https://github.com/VeckoTheGecko/xgcm", rev = "relax-pixi-build-python" } # TODO: Switch to release version after release is cut cf_xarray = "0.10.*" cftime = "1.6.*" pooch = "1.8.*" diff --git a/pyproject.toml b/pyproject.toml index 35d63ce3e..32aa744c7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,6 @@ dependencies = [ "pyarrow >=20.0.0", "uxarray >=2026.04.1", "pooch >=1.8.0", - "xgcm >=0.9.0", "cf_xarray >=0.8.6", "polars >=1.31.0", ] @@ -200,7 +199,6 @@ module = [ "cftime", "netCDF4", "pooch", - "xgcm", "uxarray", ] ignore_missing_imports = true diff --git a/src/parcels/_sgrid/__init__.py b/src/parcels/_sgrid/__init__.py index f06e03830..b0221a424 100644 --- a/src/parcels/_sgrid/__init__.py +++ b/src/parcels/_sgrid/__init__.py @@ -9,7 +9,6 @@ get_n_faces, get_n_nodes, load_mappings, - xgcm_parse_sgrid, ) __all__ = [ @@ -23,5 +22,4 @@ "get_n_faces", "get_n_nodes", "load_mappings", - "xgcm_parse_sgrid", ] diff --git a/src/parcels/_sgrid/core.py b/src/parcels/_sgrid/core.py index 05bec18ad..393e6af8d 100644 --- a/src/parcels/_sgrid/core.py +++ b/src/parcels/_sgrid/core.py @@ -38,15 +38,6 @@ class Padding(enum.Enum): BOTH = "both" -SGRID_PADDING_TO_XGCM_POSITION = { - Padding.LOW: "right", - Padding.HIGH: "left", - Padding.BOTH: "inner", - Padding.NONE: "outer", - # "center" position is not used in SGrid, in SGrid this would just be the edges/faces themselves -} - - def get_n_faces(n_nodes: int, padding: Padding) -> int: """Get number of faces along a dimension""" if padding in [Padding.LOW, Padding.HIGH]: @@ -467,31 +458,6 @@ def parse_grid_attrs(attrs: dict[str, Hashable]) -> SGrid2DMetadata | SGrid3DMet return grid -def xgcm_parse_sgrid(ds: xr.Dataset): - # Function similar to that provided in `xgcm.metadata_parsers. - # Might at some point be upstreamed to xgcm directly - grid = ds.sgrid.metadata - - if isinstance(grid, SGrid2DMetadata): - dimensions = grid.face_dimensions + (grid.vertical_dimensions or ()) - else: - assert isinstance(grid, SGrid3DMetadata) - dimensions = grid.volume_dimensions - - xgcm_coords = {} - for face_node_padding, axis in zip(dimensions, "XYZ", strict=False): - xgcm_position = SGRID_PADDING_TO_XGCM_POSITION[face_node_padding.padding] - - coords = {} - for pos, dim in [("center", face_node_padding.face), (xgcm_position, face_node_padding.node)]: - # only include dimensions in dataset (ignore dimensions in metadata that may not exist - e.g., due to `.isel`) - if dim in ds.dims: - coords[pos] = dim - xgcm_coords[axis] = coords - - return (ds, {"coords": xgcm_coords}) - - def _get_unique_names(grid: SGrid2DMetadata | SGrid3DMetadata) -> set[str]: dims = set() dims.update(set(grid.node_dimensions)) diff --git a/tests/sgrid/test_sgrid.py b/tests/sgrid/test_sgrid.py index b83bff494..f9fd3c43f 100644 --- a/tests/sgrid/test_sgrid.py +++ b/tests/sgrid/test_sgrid.py @@ -4,13 +4,12 @@ import numpy as np import pytest import xarray as xr -import xgcm from hypothesis import assume, example, given import parcels._sgrid as sgrid import parcels._strategies as pst from parcels._sgrid.accessor import get_dim_position -from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION, _get_unique_names, parse_grid_attrs +from parcels._sgrid.core import _get_unique_names, parse_grid_attrs def create_example_grid2dmetadata(with_vertical_dimensions: bool, with_node_coordinates: bool): @@ -281,46 +280,6 @@ def test_parse_grid_attrs(grid): assert parsed == grid -@example(grid2dmetadata) -@given(pst.sgrid.grid2Dmetadata()) -def test_parse_sgrid_2d(grid_metadata: sgrid.SGrid2DMetadata): - """Test the ingestion of datasets in XGCM to ensure that it matches the SGRID metadata provided""" - ds = dummy_sgrid_2d_ds(grid_metadata) - - _, xgcm_kwargs = sgrid.xgcm_parse_sgrid(ds) - grid = xgcm.Grid(ds, autoparse_metadata=False, **xgcm_kwargs) - - for obj, axis in zip(grid_metadata.face_dimensions, ["X", "Y"], strict=True): - coords = grid.axes[axis].coords - assert coords["center"] == obj.face - assert coords[SGRID_PADDING_TO_XGCM_POSITION[obj.padding]] == obj.node - - if grid_metadata.vertical_dimensions is None: - assert "Z" not in grid.axes - else: - obj = grid_metadata.vertical_dimensions[0] - coords = grid.axes["Z"].coords - assert coords["center"] == obj.face - assert coords[SGRID_PADDING_TO_XGCM_POSITION[obj.padding]] == obj.node - - -@given(pst.sgrid.grid3Dmetadata()) -@pytest.mark.xfail( - reason="Parcels doesn't have native support for SGRID 3D grids. This metadata checking is superfluous until we have such support." -) -def test_parse_sgrid_3d(grid_metadata: sgrid.SGrid3DMetadata): - """Test the ingestion of datasets in XGCM to ensure that it matches the SGRID metadata provided""" - ds = dummy_sgrid_3d_ds(grid_metadata) - - ds, xgcm_kwargs = sgrid.xgcm_parse_sgrid(ds) - grid = xgcm.Grid(ds, autoparse_metadata=False, **xgcm_kwargs) - - for obj, axis in zip(grid_metadata.volume_dimensions, ["X", "Y", "Z"], strict=True): - coords = grid.axes[axis].coords - assert coords["center"] == obj.face - assert coords[SGRID_PADDING_TO_XGCM_POSITION[obj.padding]] == obj.node - - @pytest.mark.parametrize( "grid", [ From 1aee6ffb843f5d0c2a6d2a0658ae6160674dfd3d Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 5 Aug 2026 11:01:55 +0800 Subject: [PATCH 06/13] fix: update remaining xgcm_grid references missed in initial refactor - get_cell_count_along_dim: fall back to deriving face size from node dim when the face dim is absent from the dataset - _xinterpolators.py: replace grid.xgcm_grid.axes offset checks with sgrid_metadata.face/vertical_dimensions padding comparisons - _reprs.py: replace xgcm_grid repr with sgrid_metadata repr - tests/utils.py: update get_axis_from_dim_name call to new signature --- src/parcels/_core/xgrid.py | 4 +++- src/parcels/_reprs.py | 4 ++-- src/parcels/interpolators/_xinterpolators.py | 13 +++++++------ tests/utils.py | 2 +- 4 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 345c2a211..a05b6b10f 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -19,7 +19,9 @@ def get_cell_count_along_dim(ds: xr.Dataset, fnp: sgrid.FaceNodePadding) -> int: - return ds[fnp.face].size - 1 + if fnp.face in ds.dims: + return ds.sizes[fnp.face] - 1 + return sgrid.get_n_faces(ds.sizes[fnp.node], fnp.padding) - 1 def get_time(ds: xr.Dataset, time_dim: str) -> npt.NDArray: diff --git a/src/parcels/_reprs.py b/src/parcels/_reprs.py index 5a6ae2944..56a6f53e8 100644 --- a/src/parcels/_reprs.py +++ b/src/parcels/_reprs.py @@ -79,8 +79,8 @@ def xgrid_repr(grid: Any) -> str: Parcels attributes: mesh : {grid._mesh} spatialhash : {grid._spatialhash} - xgcm Grid: -{textwrap.indent(repr(grid.xgcm_grid), 8 * " ")} + SGRID Metadata: +{textwrap.indent(repr(grid.sgrid_metadata), 8 * " ")} """ return textwrap.dedent(out).strip() diff --git a/src/parcels/interpolators/_xinterpolators.py b/src/parcels/interpolators/_xinterpolators.py index 54606c4b0..07fcebd81 100644 --- a/src/parcels/interpolators/_xinterpolators.py +++ b/src/parcels/interpolators/_xinterpolators.py @@ -9,6 +9,7 @@ from dask import is_dask_collection import parcels._core.utils.interpolation as i_u +import parcels._sgrid as sgrid import parcels._typing as ptyping from parcels.interpolators._base import ScalarInterpolator, VectorInterpolator @@ -67,12 +68,12 @@ def _get_corner_data_Agrid( def _get_offsets_dictionary(grid: XGrid) -> dict[ptyping.CfAxisSpatial, Literal[1, 0]]: offsets = {} - for axis in ["X", "Y"]: - axis_coords = grid.xgcm_grid.axes[axis].coords.keys() - offsets[axis] = 1 if "right" in axis_coords else 0 - if "Z" in grid.xgcm_grid.axes: - axis_coords = grid.xgcm_grid.axes["Z"].coords.keys() - offsets["Z"] = 1 if "right" in axis_coords else 0 + metadata = grid.sgrid_metadata + for fnp, axis in zip(metadata.face_dimensions, ["X", "Y"], strict=False): + offsets[axis] = 1 if fnp.padding == sgrid.Padding.LOW else 0 + if metadata.vertical_dimensions is not None: + fnp_z = metadata.vertical_dimensions[0] + offsets["Z"] = 1 if fnp_z.padding == sgrid.Padding.LOW else 0 else: offsets["Z"] = 0 return offsets diff --git a/tests/utils.py b/tests/utils.py index a07cf6e97..0ff8589ef 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -134,7 +134,7 @@ def assert_valid_field_data(data: xr.DataArray, grid: XGrid): assert len(data.shape) == 4, f"Field data should have 4 dimensions (time, depth, lat, lon), got dims {data.dims}" for ax_expected, dim in zip(_FIELD_DATA_ORDERING, data.dims, strict=True): - ax_actual = get_axis_from_dim_name(grid.xgcm_grid.axes, dim) + ax_actual = get_axis_from_dim_name(grid.sgrid_metadata, dim) if ax_actual is None: continue # None is ok assert ax_actual == ax_expected, f"Expected axis {ax_expected} for dimension '{dim}', got {ax_actual}" From 545742487348a34f21f86d436da3fc5a25b7ee3a Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 5 Aug 2026 11:40:38 +0800 Subject: [PATCH 07/13] Remove design documents --- .../plans/2026-08-04-remove-xgcm.md | 896 ------------------ .../specs/2026-08-04-remove-xgcm-design.md | 129 --- 2 files changed, 1025 deletions(-) delete mode 100644 docs/superpowers/plans/2026-08-04-remove-xgcm.md delete mode 100644 docs/superpowers/specs/2026-08-04-remove-xgcm-design.md diff --git a/docs/superpowers/plans/2026-08-04-remove-xgcm.md b/docs/superpowers/plans/2026-08-04-remove-xgcm.md deleted file mode 100644 index 0d52f463b..000000000 --- a/docs/superpowers/plans/2026-08-04-remove-xgcm.md +++ /dev/null @@ -1,896 +0,0 @@ -# Remove xgcm Dependency Implementation Plan - -> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. - -**Goal:** Remove xgcm as a runtime dependency by replacing the XgcmLike adapter layer with direct use of SGRID metadata throughout `xgrid.py`. - -**Architecture:** Three sequential tasks: (1) add new SGRID-native types and helpers as pure additions, (2) refactor xgrid.py to use them instead of the adapter layer, (3) delete the now-dead bridge functions and fix the two test files that called xgcm directly. - -**Tech Stack:** Python 3.11+, xarray, parcels SGRID metadata (`SGrid2DMetadata`, `FaceNodePadding`, `Padding`) - -## Global Constraints - -- Never import `xgcm` outside of `TYPE_CHECKING` blocks — and those blocks are being removed too -- Preserve all existing function names (they are imported by `model.py` and may be used elsewhere) -- All commits must include `Co-authored-by: Claude ` -- Run `pytest tests/sgrid/ tests/datasets/ -x` after each task before committing - ---- - -## File Map - -| File | Change | -| ----------------------------------- | -------------------------------------------------------------------------------------- | -| `src/parcels/_typing.py` | Add `GridPosition`; remove `XgcmAxisPosition`, `XgcmAxes`, xgcm TYPE_CHECKING import | -| `src/parcels/_sgrid/accessor.py` | Add `get_dim_position()` | -| `src/parcels/_core/xgrid.py` | Major refactor — remove adapter classes, update all functions to use `SGrid2DMetadata` | -| `src/parcels/_sgrid/core.py` | Remove `SGRID_PADDING_TO_XGCM_POSITION` and `xgcm_parse_sgrid()` | -| `src/parcels/_sgrid/__init__.py` | Remove `xgcm_parse_sgrid` from imports and `__all__` | -| `tests/sgrid/test_sgrid.py` | Remove xgcm import + two xgcm tests; add `test_get_dim_position` | -| `tests/datasets/test_structured.py` | Replace `xgcm.Grid` calls with SGRID-native assertions | -| `pyproject.toml` | Remove `"xgcm >=0.9.0"` | - ---- - -### Task 1: Add `GridPosition` type and `get_dim_position()` helper - -Pure additions — nothing is deleted. All existing code continues to work. - -**Files:** - -- Modify: `src/parcels/_typing.py` -- Modify: `src/parcels/_sgrid/accessor.py` -- Modify: `tests/sgrid/test_sgrid.py` - -**Interfaces:** - -- Produces: `GridPosition = Literal["face"] | Padding` in `_typing.py` -- Produces: `get_dim_position(grid: SGrid2DMetadata, dim: str) -> Literal["face"] | Padding` in `_sgrid/accessor.py` - -- [ ] **Step 1: Write failing tests for `get_dim_position`** - -Add to `tests/sgrid/test_sgrid.py` (after the existing imports, before existing test functions): - -```python -from parcels._sgrid.accessor import get_dim_position - - -def test_get_dim_position_face_dims(): - """Face dimensions return 'face'.""" - metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) - # face_dimensions = (FaceNodePadding("face_dimension1", "node_dimension1", Padding.LOW), ...) - assert get_dim_position(metadata, "face_dimension1") == "face" - assert get_dim_position(metadata, "face_dimension2") == "face" - - -def test_get_dim_position_node_dims(): - """Node dimensions return their Padding value.""" - metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) - assert get_dim_position(metadata, "node_dimension1") == sgrid.Padding.LOW - assert get_dim_position(metadata, "node_dimension2") == sgrid.Padding.LOW - - -def test_get_dim_position_vertical(): - """Vertical face and node dimensions are handled.""" - metadata = create_example_grid2dmetadata(with_vertical_dimensions=True, with_node_coordinates=False) - # vertical_dimensions = (FaceNodePadding("vertical_dimensions_dim1", "vertical_dimensions_dim2", Padding.LOW),) - assert get_dim_position(metadata, "vertical_dimensions_dim1") == "face" - assert get_dim_position(metadata, "vertical_dimensions_dim2") == sgrid.Padding.LOW - - -def test_get_dim_position_unknown_dim_raises(): - """Unknown dimensions raise ValueError.""" - metadata = create_example_grid2dmetadata(with_vertical_dimensions=False, with_node_coordinates=False) - with pytest.raises(ValueError, match="not a spatial SGRID dimension"): - get_dim_position(metadata, "nonexistent_dim") -``` - -- [ ] **Step 2: Run tests to confirm they fail** - -```bash -pytest tests/sgrid/test_sgrid.py::test_get_dim_position_face_dims -x -v -``` - -Expected: `ImportError` or `AttributeError` — `get_dim_position` does not exist yet. - -- [ ] **Step 3: Add `GridPosition` to `_typing.py`** - -In `src/parcels/_typing.py`, after the existing imports block, add: - -```python -from parcels._sgrid.core import Padding -``` - -Then after `XgcmAxisDirection = CfAxisSpatial | Literal["T"]`, add: - -```python -GridPosition = Literal["face"] | Padding -``` - -No existing lines are removed yet — that happens in Task 2. - -- [ ] **Step 4: Add `get_dim_position` to `accessor.py`** - -In `src/parcels/_sgrid/accessor.py`, add this function after `_get_axis_info`: - -```python -def get_dim_position(grid: SGrid2DMetadata, dim: str) -> "Literal['face'] | Padding": - """Returns 'face' if dim is a face dimension, or the Padding value if it is a node dimension. - - Replaces xgcm's position string vocabulary ('center', 'left', 'right', 'inner', 'outer') - with SGRID-native types. - """ - axis_info = _get_axis_info(grid) - if dim not in axis_info: - raise ValueError(f"Dimension {dim!r} is not a spatial SGRID dimension in this grid.") - fnp, is_node = axis_info[dim] - return fnp.padding if is_node else "face" -``` - -Add the `Literal` import to the existing `from typing import Any, Literal, cast` line in accessor.py (it already has `Literal` — confirm and leave as-is if so). - -- [ ] **Step 5: Run tests to confirm they pass** - -```bash -pytest tests/sgrid/test_sgrid.py::test_get_dim_position_face_dims tests/sgrid/test_sgrid.py::test_get_dim_position_node_dims tests/sgrid/test_sgrid.py::test_get_dim_position_vertical tests/sgrid/test_sgrid.py::test_get_dim_position_unknown_dim_raises -v -``` - -Expected: 4 PASSED. - -- [ ] **Step 6: Run full test subset to confirm no regressions** - -```bash -pytest tests/sgrid/ tests/datasets/ -x -q -``` - -Expected: all pass (same as before this task). - -- [ ] **Step 7: Commit** - -```bash -git add src/parcels/_typing.py src/parcels/_sgrid/accessor.py tests/sgrid/test_sgrid.py -git commit -m "feat: add GridPosition type and get_dim_position() SGRID helper - -Co-authored-by: Claude " -``` - ---- - -### Task 2: Refactor `xgrid.py` to use SGRID metadata directly - -Remove `XgcmLikeAxis`, `XgcmLikeGrid`, `construct_xgcm_axes_object`, and `self.xgcm_grid`. Update every function in `xgrid.py` that referenced the adapter layer to use `SGrid2DMetadata` / `FaceNodePadding` / `Padding` directly. Also clean up the now-dead type aliases from `_typing.py`. - -**Files:** - -- Modify: `src/parcels/_core/xgrid.py` -- Modify: `src/parcels/_typing.py` - -**Interfaces:** - -- Consumes: `get_dim_position(grid, dim)` from Task 1 -- Consumes: `GridPosition` from Task 1 -- Consumes: `_get_dim_to_axis_mapping(metadata)` from `_sgrid/accessor.py` (already imported) -- Produces: All public function signatures updated to take `SGrid2DMetadata` instead of `XgcmAxes`/`xgcm.Grid` - -- [ ] **Step 1: Confirm baseline tests pass before touching anything** - -```bash -pytest tests/sgrid/ tests/datasets/ -x -q -``` - -Expected: all pass. If not, do not proceed — fix the cause first. - -- [ ] **Step 2: Rewrite `xgrid.py` imports and module-level constants** - -At the top of `src/parcels/_core/xgrid.py`, make these changes: - -Remove the `TYPE_CHECKING` block entirely: - -```python -# DELETE these lines: -if TYPE_CHECKING: - import xgcm.axis -``` - -Remove the `_DEFAULT_XGCM_KWARGS` constant: - -```python -# DELETE this line: -_DEFAULT_XGCM_KWARGS: dict[str, Any] = {"padding": "fill"} -``` - -Remove the import of `SGRID_PADDING_TO_XGCM_POSITION`: - -```python -# DELETE this line: -from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION -``` - -Add `get_dim_position` to the existing accessor import: - -```python -# Change this: -from parcels._sgrid.accessor import _get_dim_to_axis_mapping -# To: -from parcels._sgrid.accessor import _get_dim_to_axis_mapping, get_dim_position -``` - -Add `FaceNodePadding` and `Padding` to the sgrid imports. The existing `import parcels._sgrid as sgrid` is kept — use `sgrid.Padding` and `sgrid.FaceNodePadding` in the code below. - -Also update the `_typing` import to add `GridPosition`: - -```python -# The existing line: -import parcels._typing as ptyping -# No change to the line itself — access GridPosition as ptyping.GridPosition -``` - -Remove `Any` from the typing imports if it's only used by `_DEFAULT_XGCM_KWARGS` (check first — if used elsewhere, keep it). - -- [ ] **Step 3: Replace standalone helper functions** - -Replace `get_cell_count_along_dim` (lines 29-33): - -```python -def get_cell_count_along_dim(ds: xr.Dataset, fnp: sgrid.FaceNodePadding) -> int: - return ds[fnp.face].size - 1 -``` - -Replace `get_time` (lines 36-37): - -```python -def get_time(ds: xr.Dataset, time_dim: str) -> npt.NDArray: - return ds[time_dim].values -``` - -Replace `_get_xgrid_axes` (lines 40-42): - -```python -def _get_xgrid_axes(metadata: sgrid.SGrid2DMetadata, ds_dims: set[str]) -> list[ptyping.XgridAxis]: - dim_to_axis = _get_dim_to_axis_mapping(metadata) - present = {axis for dim, axis in dim_to_axis.items() if dim in ds_dims} - return sorted(present, key=_XGRID_AXES_ORDERING.index) -``` - -Replace `assert_all_field_dims_have_axis` (lines 54-76). Note: `model.py` imports this by name so the name must not change: - -```python -def assert_all_field_dims_have_axis(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata) -> None: - dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} - ax_dims = [(dim_to_axis.get(str(dim)), str(dim)) for dim in da.dims] - - for ax, dim_name in ax_dims: - if ax is None: - raise ValueError( - f'Dimension "{dim_name}" has no axis attribute. ' - f'HINT: You may want to add an {{"axis": A}} to your DataSet["{dim_name}"], where A is one of "X", "Y", "Z" or "T"' - ) - - seen_axes: dict[str, str] = {} - for ax, dim_name in ax_dims: - if ax in seen_axes: - raise ValueError( - f"Two dimensions ({dim_name!r} and {seen_axes[ax]!r}) provide values in the axis direction {ax!r}. " - "This is not possible, a field cannot have two dimensions on a single axis." - ) - seen_axes[ax] = dim_name - assert len(ax_dims) <= 4, ( - "The input dataset appears to have more than 4 dimensions after conversion. Execution should never reach this point. Please file an issue sharing more about your input dataset." - ) -``` - -- [ ] **Step 4: Delete the adapter classes and construct function** - -Delete `XgcmLikeAxis` (lines 116-118), `XgcmLikeGrid` (lines 121-129), and `construct_xgcm_axes_object` (lines 132-159) entirely. - -- [ ] **Step 5: Rewrite `XGrid.__init__`** - -Replace the `__init__` body from `grid = XgcmLikeGrid(...)` onward: - -```python -def __init__(self, model_data: xr.Dataset, mesh: Literal["flat", "spherical"] | SphericalMesh): - self.sgrid_metadata = model_data.sgrid.metadata - self._ds = model_data - self._mesh = get_mesh(mesh) - self._spatialhash = None - ds = model_data - - if "lon" in ds: - ds.set_coords("lon") - if "lat" in ds: - ds.set_coords("lat") - - axes = self.axes # uses _get_xgrid_axes(self.sgrid_metadata, set(self._ds.dims)) - if len(set(axes) & {"X", "Y"}) > 0: - assert_valid_lat_lon(ds["lat"], ds["lon"], self.sgrid_metadata) - - if "Z" in axes: - assert_valid_depth(ds["depth"]) - - self._ds = ds -``` - -- [ ] **Step 6: Rewrite `XGrid` properties and methods** - -Replace the `axes` property: - -```python -@property -def axes(self) -> list[ptyping.XgridAxis]: - return _get_xgrid_axes(self.sgrid_metadata, set(self._ds.dims)) -``` - -Replace the `lon` property: - -```python -@property -def lon(self): - """ - Note - ---- - Included for compatibility with v3 codebase. May be removed in future. - TODO v4: Evaluate - """ - if "X" not in self.axes: - return np.zeros(1) - if is_dask_collection(self._ds["lon"].data): - self._ds["lon"].load() - return self._ds["lon"].values -``` - -Replace the `lat` property: - -```python -@property -def lat(self): - """ - Note - ---- - Included for compatibility with v3 codebase. May be removed in future. - TODO v4: Evaluate - """ - if "Y" not in self.axes: - return np.zeros(1) - if is_dask_collection(self._ds["lat"].data): - self._ds["lat"].load() - return self._ds["lat"].values -``` - -Replace the `depth` property: - -```python -@property -def depth(self): - """ - Note - ---- - Included for compatibility with v3 codebase. May be removed in future. - TODO v4: Evaluate - """ - if "Z" not in self.axes: - return np.zeros(1) - return self._ds["depth"].values -``` - -Replace the `_datetimes` property: - -```python -@property -def _datetimes(self): - if "time" not in self._ds.dims: - return np.zeros(1) - return get_time(self._ds, "time") -``` - -Replace `get_axis_dim`: - -```python -def get_axis_dim(self, axis: ptyping.XgridAxis) -> int: - if axis not in self.axes: - raise ValueError(f"Axis {axis!r} is not part of this grid. Available axes: {self.axes}") - - fnp_x, fnp_y = self.sgrid_metadata.face_dimensions - if axis == "X": - return get_cell_count_along_dim(self._ds, fnp_x) - if axis == "Y": - return get_cell_count_along_dim(self._ds, fnp_y) - # axis == "Z" - assert self.sgrid_metadata.vertical_dimensions is not None - return get_cell_count_along_dim(self._ds, self.sgrid_metadata.vertical_dimensions[0]) -``` - -Replace `localize`: - -```python -def localize( - self, position: dict[ptyping.XgridAxis, tuple[int, float]], dims: list[str] -) -> dict[str, tuple[int, float]]: - """ - Uses the grid context (i.e., the staggering of the grid) to convert a position relative - to the F-points in the grid to a position relative to the staggered grid the array - of interest is defined on. - - Uses dimensions of the DataArray to determine the staggered grid. - - WARNING: This API is unstable and subject to change in future versions. - - Parameters - ---------- - position : dict - A mapping of the axis to a tuple of (index, barycentric coordinate) for the - F-points in the grid. - dims : list[str] - A list of dimension names that the DataArray is defined on. This is used to determine - the staggering of the grid and which axis each dimension corresponds to. - - Returns - ------- - dict[str, tuple[int, float]] - A mapping of the dimension names to a tuple of (index, barycentric coordinate) for - the staggered grid the DataArray is defined on. - - Example - ------- - >>> position = {'X': (5, 0.51), 'Y': ( - 10, 0.25), 'Z': (3, 0.75)} - >>> dims = ['time', 'depth', 'YC', 'XC'] - >>> grid.localize(position, dims) - {'depth': (3, 0.75), 'YC': (9, 0.75), 'XC': (5, 0.01)} - """ - dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) | {"time": "T"} - axis_to_var = {dim_to_axis[dim]: dim for dim in dims if dim in dim_to_axis} - var_positions = { - axis: get_dim_position(self.sgrid_metadata, dim) - for axis, dim in axis_to_var.items() - if axis != "T" - } - return { - axis_to_var[axis]: _convert_center_pos_to_fpoint( - index=index, - bcoord=bcoord, - position=var_positions[axis], - f_point_position=self._fpoint_info[axis], - ) - for axis, (index, bcoord) in position.items() - } -``` - -Replace `_fpoint_info`: - -```python -@cached_property -def _fpoint_info(self) -> dict[ptyping.XgridAxis, sgrid.Padding]: - """Returns a mapping of the spatial axes in the Grid to their Padding values (node positions).""" - metadata = self.sgrid_metadata - fnp_x, fnp_y = metadata.face_dimensions - result: dict[ptyping.XgridAxis, sgrid.Padding] = {} - axes = self.axes - if "X" in axes: - result["X"] = fnp_x.padding - if "Y" in axes: - result["Y"] = fnp_y.padding - if "Z" in axes and metadata.vertical_dimensions: - result["Z"] = metadata.vertical_dimensions[0].padding - return result -``` - -Replace `get_axis_dim_mapping`: - -```python -def get_axis_dim_mapping(self, dims: Sequence[Hashable]) -> dict[ptyping.XgridAxis, str]: - """ - Maps xarray dimension names to their corresponding axis (X, Y, Z). - - WARNING: This API is unstable and subject to change in future versions. - - Parameters - ---------- - dims : Sequence[Hashable] - Sequence of xarray dimension names - - Returns - ------- - dict[_XGRID_AXES, str] - Dictionary mapping axes (X, Y, Z) to their corresponding dimension names - - Examples - -------- - >>> grid.get_axis_dim_mapping(['time', 'lat', 'lon']) - {'Y': 'lat', 'X': 'lon'} - - Notes - ----- - Only returns mappings for spatial axes (X, Y, Z) that are present in the grid. - """ - dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) - result = {} - for dim in dims: - axis = dim_to_axis.get(str(dim)) - if axis in self.axes: - result[cast(ptyping.XgridAxis, axis)] = str(dim) - return result -``` - -- [ ] **Step 7: Rewrite module-level functions below `XGrid`** - -Replace `get_axis_from_dim_name` (lines 459-464): - -```python -def get_axis_from_dim_name(metadata: sgrid.SGrid2DMetadata, dim: Hashable) -> ptyping.XgcmAxisDirection | None: - """For a given dimension name in a grid, returns the direction axis it is on.""" - dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} - return dim_to_axis.get(str(dim)) -``` - -Replace `get_xgcm_position_from_dim_name` (lines 467-473). Keep the name as-is since it may be imported externally: - -```python -def get_xgcm_position_from_dim_name(metadata: sgrid.SGrid2DMetadata, dim: str) -> ptyping.GridPosition | None: - """For a given dimension, returns the GridPosition of the variable in the grid.""" - try: - return get_dim_position(metadata, dim) - except ValueError: - return None -``` - -Replace `assert_all_dimensions_correspond_with_axis` (lines 477-484): - -```python -def assert_all_dimensions_correspond_with_axis(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata) -> None: - dim_to_axis = _get_dim_to_axis_mapping(metadata) - for dim in da.dims: - if dim not in dim_to_axis: - raise ValueError( - f"Dimension {dim!r} for DataArray {da.name!r} with dims {da.dims} is not associated with a direction on the provided grid." - ) -``` - -Replace `assert_valid_field_array` (lines 487-509): - -```python -def assert_valid_field_array(da: xr.DataArray, metadata: sgrid.SGrid2DMetadata): - """ - Asserts that for a data array: - - All dimensions are associated with a direction on the grid - - These directions are T, Z, Y, X and the array is ordered as T, Z, Y, X - """ - dim_to_axis = _get_dim_to_axis_mapping(metadata) | {"time": "T"} - - for dim in da.dims: - if dim not in dim_to_axis: - raise ValueError( - f"Dimension {dim!r} for DataArray {da.name!r} with dims {da.dims} is not associated with a direction on the provided grid." - ) - - dim_to_axis_for_da = {dim: dim_to_axis[dim] for dim in da.dims} - dim_to_axis_for_da = cast(dict[Hashable, ptyping.XgcmAxisDirection], dim_to_axis_for_da) - - if set(dim_to_axis_for_da.values()) != {"T", "Z", "Y", "X"}: - raise ValueError( - f"DataArray {da.name!r} with dims {da.dims} has directions {tuple(dim_to_axis_for_da.values())}." - "Expected directions of 'T', 'Z', 'Y', and 'X'." - ) - - if list(dim_to_axis_for_da.values()) != ["T", "Z", "Y", "X"]: - raise ValueError( - f"Dimension order for array {da.name!r} is not valid. Got {tuple(dim_to_axis_for_da.keys())} with associated directions of {tuple(dim_to_axis_for_da.values())}. Expected directions of ('T', 'Z', 'Y', 'X'). Transpose your array accordingly." - ) -``` - -Replace `assert_valid_lat_lon` (lines 512-579): - -```python -def assert_valid_lat_lon(da_lat, da_lon, metadata: sgrid.SGrid2DMetadata): - """ - Asserts that the provided longitude and latitude DataArrays are defined appropriately - on the F points to match the internal representation in Parcels. - - - Longitude and latitude must be 1D or 2D (both must have the same dimensionality) - - Both are defined on the node points (i.e., not the face/center) - - If 1D: - - Longitude is associated with the X axis - - Latitude is associated with the Y axis - - If 2D: - - Lon and lat are defined on the same dimensions - - Lon and lat are transposed such they're Y, X - """ - assert_all_dimensions_correspond_with_axis(da_lon, metadata) - assert_all_dimensions_correspond_with_axis(da_lat, metadata) - - for dim in da_lon.dims: - if get_dim_position(metadata, dim) == "face": - raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is defined on the center of the grid, but must be defined on the F points." - ) - for dim in da_lat.dims: - if get_dim_position(metadata, dim) == "face": - raise ValueError( - f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is defined on the center of the grid, but must be defined on the F points." - ) - - if da_lon.ndim != da_lat.ndim: - raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} have different dimensionalities." - ) - if da_lon.ndim not in (1, 2): - raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be 1D or 2D." - ) - - dim_to_axis = _get_dim_to_axis_mapping(metadata) - - if da_lon.ndim == 1: - if dim_to_axis.get(da_lon.dims[0]) != "X": - raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is not associated with the X axis." - ) - if dim_to_axis.get(da_lat.dims[0]) != "Y": - raise ValueError( - f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is not associated with the Y axis." - ) - - if not np.all(np.diff(da_lon.values) > 0): - raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} must be strictly increasing." - ) - if not np.all(np.diff(da_lat.values) > 0): - raise ValueError(f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be strictly increasing.") - - if da_lon.ndim == 2: - if da_lon.dims != da_lat.dims: - raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be defined on the same dimensions." - ) - - lon_axes = [dim_to_axis.get(dim) for dim in da_lon.dims] - if lon_axes != ["Y", "X"]: - raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} and Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} must be defined on the X and Y axes and transposed to have dimensions in order of Y, X." - ) -``` - -Replace `_convert_center_pos_to_fpoint` (lines 590-616): - -```python -def _convert_center_pos_to_fpoint( - *, - index: int, - bcoord: float, - position: ptyping.GridPosition, - f_point_position: sgrid.Padding, -) -> tuple[int, float]: - """Converts a physical position relative to the cell edges defined in the grid to be relative to the center point. - - This is used to "localize" a position to be relative to the staggered grid at which the field is defined, so that - it can be easily interpolated. - - This also handles different model input cell edges and centers are staggered in different directions (e.g., with NEMO and MITgcm). - """ - if position != "face": # Data is already defined on the F points - return index, bcoord - - bcoord = bcoord - 0.5 - if bcoord < 0: - bcoord += 1.0 - index -= 1 - - # Correct relative to the f-point position - # Padding.BOTH was "inner", Padding.LOW was "right" in xgcm vocabulary - if f_point_position in (sgrid.Padding.BOTH, sgrid.Padding.LOW): - index += 1 - - return index, bcoord -``` - -- [ ] **Step 8: Clean up `_typing.py`** - -Remove these lines from `src/parcels/_typing.py`: - -```python -# DELETE the TYPE_CHECKING block for xgcm: -if TYPE_CHECKING: - import xgcm - -# DELETE these two type aliases: -XgcmAxisPosition = Literal["center", "left", "right", "inner", "outer"] -XgcmAxes = Mapping[XgcmAxisDirection, "xgcm.Axis"] -``` - -Also remove `Mapping` from the `collections.abc` import if it is now unused (check by searching for other uses of `Mapping` in the file). - -- [ ] **Step 9: Run tests** - -```bash -pytest tests/sgrid/ tests/datasets/ -x -q -``` - -Expected: all pass. If any fail, fix before proceeding. - -- [ ] **Step 10: Commit** - -```bash -git add src/parcels/_core/xgrid.py src/parcels/_typing.py -git commit -m "refactor: replace XgcmLike adapter layer with direct SGRID metadata usage - -Co-authored-by: Claude " -``` - ---- - -### Task 3: Remove bridge functions, fix tests, remove xgcm dep - -Delete the now-dead bridge code in `_sgrid/core.py`, fix the two test files, and drop xgcm from `pyproject.toml`. - -**Files:** - -- Modify: `src/parcels/_sgrid/core.py` -- Modify: `src/parcels/_sgrid/__init__.py` -- Modify: `tests/sgrid/test_sgrid.py` -- Modify: `tests/datasets/test_structured.py` -- Modify: `pyproject.toml` - -**Interfaces:** - -- Consumes: `get_dim_position` from Task 1 - -- [ ] **Step 1: Remove `SGRID_PADDING_TO_XGCM_POSITION` and `xgcm_parse_sgrid` from `_sgrid/core.py`** - -Delete lines 41-47 (the `SGRID_PADDING_TO_XGCM_POSITION` dict): - -```python -# DELETE: -SGRID_PADDING_TO_XGCM_POSITION = { - Padding.LOW: "right", - Padding.HIGH: "left", - Padding.BOTH: "inner", - Padding.NONE: "outer", - # "center" position is not used in SGrid, in SGrid this would just be the edges/faces themselves -} -``` - -Delete lines 470-492 (the `xgcm_parse_sgrid` function): - -```python -# DELETE: -def xgcm_parse_sgrid(ds: xr.Dataset): - # Function similar to that provided in `xgcm.metadata_parsers. - # Might at some point be upstreamed to xgcm directly - grid = ds.sgrid.metadata - ... - return (ds, {"coords": xgcm_coords}) -``` - -- [ ] **Step 2: Update `_sgrid/__init__.py`** - -Remove `xgcm_parse_sgrid` from the import and `__all__`: - -```python -# Change from: -from .core import ( - FaceNodePadding, - Padding, - SGrid2DMetadata, - SGrid3DMetadata, - _attach_sgrid_metadata, - dump_mappings, - get_n_faces, - get_n_nodes, - load_mappings, - xgcm_parse_sgrid, -) - -__all__ = [ - "FaceNodePadding", - "Padding", - "SGrid2DMetadata", - "SGrid3DMetadata", - "SgridAccessor", - "_attach_sgrid_metadata", - "dump_mappings", - "get_n_faces", - "get_n_nodes", - "load_mappings", - "xgcm_parse_sgrid", -] - -# To: -from .core import ( - FaceNodePadding, - Padding, - SGrid2DMetadata, - SGrid3DMetadata, - _attach_sgrid_metadata, - dump_mappings, - get_n_faces, - get_n_nodes, - load_mappings, -) - -__all__ = [ - "FaceNodePadding", - "Padding", - "SGrid2DMetadata", - "SGrid3DMetadata", - "SgridAccessor", - "_attach_sgrid_metadata", - "dump_mappings", - "get_n_faces", - "get_n_nodes", - "load_mappings", -] -``` - -- [ ] **Step 3: Fix `tests/sgrid/test_sgrid.py`** - -Remove the `import xgcm` line (line 7). - -Remove `SGRID_PADDING_TO_XGCM_POSITION` from the import on line 12: - -```python -# Change from: -from parcels._sgrid.core import SGRID_PADDING_TO_XGCM_POSITION, _get_unique_names, parse_grid_attrs - -# To: -from parcels._sgrid.core import _get_unique_names, parse_grid_attrs -``` - -Delete the entire `test_parse_sgrid_2d` function (lines 253-273) and `test_parse_sgrid_3d` function (lines 276-290) — these test the now-deleted `xgcm_parse_sgrid` bridge function. The SGRID parsing itself is still exercised by the existing `test_Grid2DMetadata_roundtrip`, `test_parse_grid_attrs`, and other tests. - -- [ ] **Step 4: Fix `tests/datasets/test_structured.py`** - -Replace the entire file content: - -```python -import parcels._sgrid as sgrid -from parcels._datasets.structured.generic import datasets -from parcels._sgrid.accessor import get_dim_position - - -def test_left_indexed_dataset(): - """Checks that 'ds_2d_left' uses HIGH padding (MITgcm / left-indexed style).""" - ds = datasets["ds_2d_left"] - metadata = ds.sgrid.metadata - for fnp in metadata.face_dimensions: - assert get_dim_position(metadata, fnp.face) == "face" - assert get_dim_position(metadata, fnp.node) == sgrid.Padding.HIGH - - -def test_right_indexed_dataset(): - """Checks that 'ds_2d_right' uses LOW padding (NEMO / right-indexed style).""" - ds = datasets["ds_2d_right"] - metadata = ds.sgrid.metadata - for fnp in metadata.face_dimensions: - assert get_dim_position(metadata, fnp.face) == "face" - assert get_dim_position(metadata, fnp.node) == sgrid.Padding.LOW -``` - -- [ ] **Step 5: Remove xgcm from `pyproject.toml`** - -Delete the line `"xgcm >=0.9.0",` from the `dependencies` list in `pyproject.toml`. - -- [ ] **Step 6: Run the full test subset** - -```bash -pytest tests/sgrid/ tests/datasets/ -x -q -``` - -Expected: all pass. The two deleted tests (`test_parse_sgrid_2d`, `test_parse_sgrid_3d`) are gone; four new tests from Task 1 and two rewritten tests from `test_structured.py` all pass. - -- [ ] **Step 7: Verify xgcm is no longer imported anywhere in src/** - -```bash -grep -r "import xgcm" src/parcels --include="*.py" -``` - -Expected: no output. - -- [ ] **Step 8: Commit** - -```bash -git add src/parcels/_sgrid/core.py src/parcels/_sgrid/__init__.py \ - tests/sgrid/test_sgrid.py tests/datasets/test_structured.py \ - pyproject.toml -git commit -m "feat: remove xgcm dependency — use SGRID metadata natively throughout - -Co-authored-by: Claude " -``` diff --git a/docs/superpowers/specs/2026-08-04-remove-xgcm-design.md b/docs/superpowers/specs/2026-08-04-remove-xgcm-design.md deleted file mode 100644 index 310b53961..000000000 --- a/docs/superpowers/specs/2026-08-04-remove-xgcm-design.md +++ /dev/null @@ -1,129 +0,0 @@ -# Remove xgcm Dependency — Design Spec - -**Date:** 2026-08-04 -**Branch:** remove-xgcm -**Approach:** Option B — full SGRID-native refactor (remove XgcmLike* adapter layer) - -## Motivation - -All grid topology information previously sourced from xgcm (axis directions, staggering positions, face/node relationships) is now fully expressed in SGRID metadata attached to the dataset. The XgcmLike* adapter layer was always a temporary bridge. Removing xgcm simplifies the dependency graph, removes COMODO metadata coupling, and lets the codebase work entirely with the standardised SGRID model. - ---- - -## Section 1: New Position Vocabulary - -Replace xgcm position strings (`"center"`, `"left"`, `"right"`, `"inner"`, `"outer"`) with a new `GridPosition` type in `_typing.py`: - -```python -from parcels._sgrid.core import Padding - -GridPosition = Literal["face"] | Padding -``` - -| Old (xgcm) | New (SGRID) | -| ---------- | -------------- | -| `"center"` | `"face"` | -| `"right"` | `Padding.LOW` | -| `"left"` | `Padding.HIGH` | -| `"inner"` | `Padding.BOTH` | -| `"outer"` | `Padding.NONE` | - -**Removed from `_typing.py`:** - -- `XgcmAxisPosition` -- `XgcmAxes` -- `import xgcm` (TYPE_CHECKING block) - -**Removed from `_sgrid/core.py`:** - -- `SGRID_PADDING_TO_XGCM_POSITION` dict -- `xgcm_parse_sgrid()` function - -**Added to `_sgrid/accessor.py`:** - -```python -def get_dim_position(grid: SGrid2DMetadata, dim: str) -> GridPosition: - """Returns 'face' or the Padding value for a given dimension.""" -``` - -This replaces `get_xgcm_position_from_dim_name`. Uses the existing `_get_axis_info` helper internally. - ---- - -## Section 2: `xgrid.py` Changes - -### Removed entirely - -- `XgcmLikeAxis` dataclass -- `XgcmLikeGrid` class -- `construct_xgcm_axes_object` function -- `self.xgcm_grid` attribute on `XGrid` -- `_DEFAULT_XGCM_KWARGS` -- `import xgcm` and `import xgcm.axis` (TYPE_CHECKING blocks) -- Import of `SGRID_PADDING_TO_XGCM_POSITION` - -`self.sgrid_metadata` (already stored on `XGrid`) becomes the sole source of grid topology truth. - -### Changed function signatures - -| Function | Old signature | New signature | -| -------------------------------------------- | ----------------------------------------------------------------------------- | ------------------------------------------------------- | -| `get_cell_count_along_dim` | `(ds, axis: xgcm.axis.Axis)` | `(ds, fnp: FaceNodePadding)` — uses `ds[fnp.face].size` | -| `get_time` | `(ds, axis: xgcm.axis.Axis)` | `(ds, time_dim: str)` — uses `ds[time_dim].values` | -| `_get_xgrid_axes` | `(grid: xgcm.Grid)` | `(metadata: SGrid2DMetadata)` | -| `assert_all_field_dims_have_axis` | `(da, xgcm_grid: xgcm.Grid)` | `(da, metadata: SGrid2DMetadata)` | -| `assert_valid_lat_lon` | `(da_lat, da_lon, axes: XgcmAxes)` | `(da_lat, da_lon, metadata: SGrid2DMetadata)` | -| `assert_all_dimensions_correspond_with_axis` | `(da, axes: XgcmAxes)` | `(da, metadata: SGrid2DMetadata)` | -| `assert_valid_field_array` | `(da, axes: XgcmAxes)` | `(da, metadata: SGrid2DMetadata)` | -| `_convert_center_pos_to_fpoint` | `(xgcm_position: XgcmAxisPosition, f_points_xgcm_position: XgcmAxisPosition)` | `(position: GridPosition, f_point_position: Padding)` | - -### Internal logic changes - -- `get_axis_from_dim_name(axes, dim)` → replaced by `_get_dim_to_axis_mapping(metadata).get(dim)` from accessor -- `get_xgcm_position_from_dim_name(axes, dim)` → replaced by `get_dim_position(metadata, dim)` from accessor -- `_fpoint_info` → returns `dict[XgridAxis, Padding]` (was `dict[XgridAxis, str]` with xgcm strings) -- `localize()` → uses `self.sgrid_metadata` instead of `self.xgcm_grid` -- `get_axis_dim_mapping()` → uses `_get_dim_to_axis_mapping(self.sgrid_metadata)` directly -- `_convert_center_pos_to_fpoint`: `"center"` branch becomes `"face"` check; `"inner"/"right"` checks become `Padding.BOTH / Padding.LOW` -- `XGrid.lon/lat/depth/_datetimes` — currently check `self.xgcm_grid.axes["X"/"Y"/"Z"/"T"]` to determine axis presence; replaced by checking `_get_dim_to_axis_mapping(self.sgrid_metadata)` for spatial axes, and `"time" in self._ds.dims` for the time axis - ---- - -## Section 3: Test Changes - -### `tests/datasets/test_structured.py` - -Replace `xgcm.Grid` calls with SGRID-native checks: - -```python -# Old -grid = xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS) -for _axis_name, axis in grid.axes.items(): - for pos, _dim_name in axis.coords.items(): - assert pos in ["left", "center"] - -# New -metadata = ds.sgrid.metadata -for fnp in metadata.face_dimensions: - assert get_dim_position(metadata, fnp.face) == "face" - assert get_dim_position(metadata, fnp.node) in (Padding.HIGH, Padding.LOW, Padding.BOTH, Padding.NONE) -``` - -The specific padding assertion per test depends on the dataset fixture (`ds_2d_left` vs `ds_2d_right`). - -### `tests/sgrid/test_sgrid.py` - -- Delete the two tests at lines 259–285 that call `xgcm_parse_sgrid()` and `xgcm.Grid(...)` — they were testing the now-removed bridge function -- Remove `import xgcm` and `SGRID_PADDING_TO_XGCM_POSITION` from imports - -### `pyproject.toml` - -Remove `"xgcm >=0.9.0"` from the dependencies list. - ---- - -## Out of Scope - -- Changing interpolation logic or search algorithms in `xgrid.py` -- Renaming `XgcmAxisDirection` / `CfAxis` type aliases (they don't reference xgcm at runtime) -- Any changes to SGRID parsing or accessor logic beyond adding `get_dim_position` From cff1cf0f044b6aa3a4a10b869eb1ec3fe07b0a06 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:42:10 +0800 Subject: [PATCH 08/13] Remaining XGCM cleanup --- .github/ci/policy.yaml | 1 - .github/ci/recipe.yaml | 1 - docs/user_guide/examples/explanation_grids.md | 1 - .../getting_started/explanation_concepts.md | 2 +- src/parcels/_core/xgrid.py | 2 +- src/parcels/_sgrid/accessor.py | 6 +----- tests/test_xgrid.py | 13 +------------ 7 files changed, 4 insertions(+), 22 deletions(-) diff --git a/.github/ci/policy.yaml b/.github/ci/policy.yaml index 2d7b6f2e5..0b5fe6371 100644 --- a/.github/ci/policy.yaml +++ b/.github/ci/policy.yaml @@ -19,7 +19,6 @@ policy: - pytest-cov - nbval - uxarray # undergoes active development - - xgcm # undergoes active development # these packages don't fail the CI, but will be printed in the report ignored_violations: [] diff --git a/.github/ci/recipe.yaml b/.github/ci/recipe.yaml index fe9f08dfa..8d270376b 100644 --- a/.github/ci/recipe.yaml +++ b/.github/ci/recipe.yaml @@ -39,7 +39,6 @@ requirements: - pandas >=2.2 - pyarrow >=20.0.0 - cf_xarray >=0.8.6 - - xgcm >=0.9.0 - zarr >=3 - uxarray>=2026.04.1 - pyogrio # needed for geopandas (uxarray -> geoviews -> geopandas -> pyogrio, but for some reason conda doesn't pick it up automatically) diff --git a/docs/user_guide/examples/explanation_grids.md b/docs/user_guide/examples/explanation_grids.md index 61484bdbf..6542f26fb 100644 --- a/docs/user_guide/examples/explanation_grids.md +++ b/docs/user_guide/examples/explanation_grids.md @@ -5,6 +5,5 @@ Parcels `Field` objects exist on a (structured) `parcels.XGrid` or (unstructured ```{note} The contents for this page are still under development in v4. TODO -- link to xgcm.Grid documentation - adapt from v3 grid indexing tutorial (../examples_v3/documentation_indexing.ipynb) ``` diff --git a/docs/user_guide/getting_started/explanation_concepts.md b/docs/user_guide/getting_started/explanation_concepts.md index 107e583a7..a83b229ca 100644 --- a/docs/user_guide/getting_started/explanation_concepts.md +++ b/docs/user_guide/getting_started/explanation_concepts.md @@ -52,7 +52,7 @@ fieldset += parcels.FieldSet.from_sgrid_conventions(ds_fset, vector_fields={"UVs ### Grid -Each `parcels.Field` is defined on a grid. With Parcels, we can simulate particles in fields on both structured (**`parcels.XGrid`**) and unstructured (**`parcels.UxGrid`**) grids. The grid is defined by the coordinates of grid cell nodes, edges, and faces. `parcels.XGrid` objects are based on [`xgcm.Grid`](https://xgcm.readthedocs.io/en/latest/grids.html), while `parcels.UxGrid` objects are based on [`uxarray.Grid`](https://uxarray.readthedocs.io/en/stable/generated/uxarray.Grid.html#uxarray.Grid) objects. +Each `parcels.Field` is defined on a grid. With Parcels, we can simulate particles in fields on both structured (**`parcels.XGrid`**) and unstructured (**`parcels.UxGrid`**) grids. The grid is defined by the coordinates of grid cell nodes, edges, and faces. `parcels.XGrid` objects are based on Xarray Datasets with attached SGRID metadata, while `parcels.UxGrid` objects are based on [`uxarray.Grid`](https://uxarray.readthedocs.io/en/stable/generated/uxarray.Grid.html#uxarray.Grid) objects. ```{admonition} 📖 Read more about grids :class: seealso diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index a05b6b10f..e7e81f777 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -70,7 +70,7 @@ def assert_all_field_dims_have_axis(da: xr.DataArray, metadata: sgrid.SGrid2DMet def _transpose_xfield_data_to_tzyx(da: xr.DataArray, sgrid_metadata: sgrid.SGrid2DMetadata) -> xr.DataArray: """ - Transpose a DataArray of any shape into a 4D array of order TZYX. Uses xgcm to determine + Transpose a DataArray of any shape into a 4D array of order TZYX. Uses SGRID metadata to determine the axes, and inserts mock dimensions of size 1 for any axes not present in the DataArray. """ dim_to_axis = _get_dim_to_axis_mapping(sgrid_metadata) | {"time": "T"} diff --git a/src/parcels/_sgrid/accessor.py b/src/parcels/_sgrid/accessor.py index 62fede102..04d050d74 100644 --- a/src/parcels/_sgrid/accessor.py +++ b/src/parcels/_sgrid/accessor.py @@ -149,11 +149,7 @@ def _get_axis_info(grid: SGrid2DMetadata) -> dict[str, tuple[FaceNodePadding, bo def get_dim_position(grid: SGrid2DMetadata, dim: str) -> "Literal['face'] | Padding": - """Returns 'face' if dim is a face dimension, or the Padding value if it is a node dimension. - - Replaces xgcm's position string vocabulary ('center', 'left', 'right', 'inner', 'outer') - with SGRID-native types. - """ + """Returns 'face' if dim is a face dimension, or the SGRID Padding value if it is a node dimension.""" axis_info = _get_axis_info(grid) if dim not in axis_info: raise ValueError(f"Dimension {dim!r} is not a spatial SGRID dimension in this grid.") diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py index d957ab0b9..c22437281 100644 --- a/tests/test_xgrid.py +++ b/tests/test_xgrid.py @@ -6,7 +6,7 @@ import xarray as xr from numpy.testing import assert_allclose -from parcels import Field, FieldSet +from parcels import FieldSet from parcels._core.index_search import ( LEFT_OUT_OF_BOUNDS, RIGHT_OUT_OF_BOUNDS, @@ -18,7 +18,6 @@ _transpose_xfield_data_to_tzyx, ) from parcels._datasets.structured.generic import X, Y, Z, datasets, datasets_sgrid -from parcels.interpolators import XLinear from tests import utils GridTestCase = namedtuple("GridTestCase", ["ds", "attr", "expected"]) @@ -148,16 +147,6 @@ def test_invalid_depth(): XGrid.from_dataset(ds, mesh="flat") -@pytest.mark.skip( - "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646" -) # TODO: axis checking no longer relies on these axis attributes being set (since we inspect the sgrid metadata directly) - I think this might be able to be removed entirely since sgrid metadata have quite informative error messaging. For planned future PR that deals with xgcm related cleanup -def test_dim_without_axis(): - ds = xr.Dataset({"z1d": (["depth"], [0])}, coords={"depth": [0]}) - grid = XGrid.from_dataset(ds, mesh="flat") - with pytest.raises(ValueError, match='Dimension "depth" has no axis attribute*'): - Field("z1d", ds["z1d"], grid, XLinear) - - @pytest.mark.skip( "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646" ) # TODO: I think we can just rely on the SGRID metadata for this (which already has robust error messaging). How should discrepencies between SGRID and axis attr be handled? From 019a8c0672952671675fe666a1dc0a4088290ba4 Mon Sep 17 00:00:00 2001 From: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 6 Aug 2026 08:32:20 +0800 Subject: [PATCH 09/13] Update src/parcels/_core/xgrid.py Co-authored-by: Erik van Sebille --- src/parcels/_core/xgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index e7e81f777..92d27fa3f 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -478,7 +478,7 @@ def assert_valid_lat_lon(da_lat, da_lon, metadata: sgrid.SGrid2DMetadata): for dim in da_lon.dims: if get_dim_position(metadata, dim) == "face": raise ValueError( - f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is defined on the center of the grid, but must be defined on the F points." + f"Longitude DataArray {da_lon.name!r} with dims {da_lon.dims} is defined on the faces of the grid, but must be defined on the F nodes." ) for dim in da_lat.dims: if get_dim_position(metadata, dim) == "face": From d695ee961b187e3c6906fbae5cfa3e053e3d2916 Mon Sep 17 00:00:00 2001 From: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 6 Aug 2026 08:33:56 +0800 Subject: [PATCH 10/13] Update src/parcels/_core/xgrid.py Co-authored-by: Erik van Sebille --- src/parcels/_core/xgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 92d27fa3f..21b94787b 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -483,7 +483,7 @@ def assert_valid_lat_lon(da_lat, da_lon, metadata: sgrid.SGrid2DMetadata): for dim in da_lat.dims: if get_dim_position(metadata, dim) == "face": raise ValueError( - f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is defined on the center of the grid, but must be defined on the F points." + f"Latitude DataArray {da_lat.name!r} with dims {da_lat.dims} is defined on the faces of the grid, but must be defined on the F nodes." ) if da_lon.ndim != da_lat.ndim: From 21ec4024007aa9adea940af5cf6c7465ac08650c Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 6 Aug 2026 08:43:43 +0800 Subject: [PATCH 11/13] Review feedback --- src/parcels/_core/xgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 21b94787b..706b01f49 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -265,7 +265,7 @@ def localize( >>> grid.localize(position, dims) {'depth': (3, 0.75), 'YC': (9, 0.75), 'XC': (5, 0.01)} """ - dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) | {"time": "T"} + dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) axis_to_var = {dim_to_axis[dim]: dim for dim in dims if dim in dim_to_axis} var_positions = { axis: get_dim_position(self.sgrid_metadata, dim) for axis, dim in axis_to_var.items() if axis != "T" From 66b0fde946e046043d43e82c742445c8e311bcfd Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 6 Aug 2026 08:46:56 +0800 Subject: [PATCH 12/13] Review feedback --- src/parcels/_core/xgrid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 706b01f49..4d330a61e 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -542,14 +542,14 @@ def _convert_center_pos_to_fpoint( position: ptyping.GridPosition, f_point_position: sgrid.Padding, ) -> tuple[int, float]: - """Converts a physical position relative to the cell edges defined in the grid to be relative to the center point. + """Converts a physical position relative to the cell edges defined in the grid to be relative to the face center. This is used to "localize" a position to be relative to the staggered grid at which the field is defined, so that it can be easily interpolated. This also handles different model input cell edges and centers are staggered in different directions (e.g., with NEMO and MITgcm). """ - if position != "face": # Data is already defined on the F points + if position != "face": # Data is already defined on the F nodes return index, bcoord bcoord = bcoord - 0.5 From ff1c653c6615130e2b58b2d68d9ce4c53b067b35 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Thu, 6 Aug 2026 14:03:55 +0800 Subject: [PATCH 13/13] Review feedback --- src/parcels/_core/xgrid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 4d330a61e..70603d914 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -235,7 +235,7 @@ def localize( ) -> dict[str, tuple[int, float]]: """ Uses the grid context (i.e., the staggering of the grid) to convert a position relative - to the F-points in the grid to a position relative to the staggered grid the array + to the f-points in the grid to a position relative to the staggered grid the array of interest is defined on. Uses dimensions of the DataArray to determine the staggered grid. @@ -246,7 +246,7 @@ def localize( ---------- position : dict A mapping of the axis to a tuple of (index, barycentric coordinate) for the - F-points in the grid. + f-points in the grid. dims : list[str] A list of dimension names that the DataArray is defined on. This is used to determine the staggering of the grid and which axis each dimension corresponds to.