Add error message for >4D field input

src/parcels/_core/xgrid.py

@@ -52,6 +52,21 @@ def assert_all_field_dims_have_axis(da: xr.DataArray, xgcm_grid: xgcm.Grid) -> N
                 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"'
             )
 
+    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:
+            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."
+    )
+    return
+
 
 def _transpose_xfield_data_to_tzyx(da: xr.DataArray, xgcm_grid: xgcm.Grid) -> xr.DataArray:
     """

tests/test_xgrid.py

@@ -6,7 +6,7 @@
 import xarray as xr
 from numpy.testing import assert_allclose
 
-from parcels import Field
+from parcels import Field, FieldSet
 from parcels._core.index_search import (
     LEFT_OUT_OF_BOUNDS,
     RIGHT_OUT_OF_BOUNDS,
@@ -17,7 +17,7 @@
     XGrid,
     _transpose_xfield_data_to_tzyx,
 )
-from parcels._datasets.structured.generic import X, Y, Z, datasets
+from parcels._datasets.structured.generic import X, Y, Z, datasets, datasets_sgrid
 from parcels.interpolators import XLinear
 from tests import utils
 
@@ -151,6 +151,29 @@ def test_dim_without_axis():
         Field("z1d", ds["z1d"], grid, XLinear)
 
 
+def test_dim_with_duplicate_axis():
+    ds = datasets_sgrid["ds_2d_padded_low"].copy()
+
+    # Add an extra Z axis
+    ds = ds[["data_g", "grid"]]
+    ds["data_g"] = ds["data_g"].expand_dims("vertical_dimensions_dim2", 1)
+
+    z = ds["vertical_dimensions_dim2"]
+    z.attrs.update({"axis": "Z", "c_grid_axis_shift": 0})
+    ds["vertical_dimensions_dim2"] = z
+
+    # TODO: Clean up this attribute setting (really, this should be on the source datasets)
+    lon = ds["lon"]
+    lat = ds["lat"]
+    lon.attrs.update({"units": "metres"})
+    lat.attrs.update({"units": "metres"})
+    ds["lon"] = lon
+    ds["lat"] = lat
+
+    with pytest.raises(ValueError, match="Two dimensions .*provide values in the axis direction 'Z'."):
+        FieldSet.from_sgrid_conventions(ds)
+
+
 def test_vertical1D_field():
     nz = 11
     ds = xr.Dataset(