Error when number of vertical dimensions is not 1 or 4 by erikvansebille · Pull Request #3013 · Parcels-code/Parcels

erikvansebille · 2026-10-26T09:50:53Z

Description

Following the comment at #2954 (comment), this PR moves the check on number of vertical dimensions from index_search to the creation of the FieldSet.

Comment by @VeckoTheGecko

I think your proposal of setting the depth as a coordinate is best (in the user script, but ultimately in our ingestion code of model.from_sgrid_conventions. I don't think putting checks in _search_1d_array makes sense since its checking something that we never expect to happen in our internal model (hence is wasted compute on every index search, or confusing to future devs).

Checklist

Improves on Error when calling _search_1d_array with a multidimensional array #2956
Tests added
This PR targets the correct branch (main for normal development, v3-support for v3 support)

…eldset creation

VeckoTheGecko

Minor comment about test location and name

wyatt-fluidnumerics

This looks good, we probably want to add a similar check in from_ugrid_convention, but that can be a separate PR.

VeckoTheGecko · 2026-10-26T14:15:17Z

+        assert_1D_or4D_vertical_dimensions(ds)
+


This looks good, we probably want to add a similar check in from_ugrid_convention, but that can be a separate PR.

It's here in this PR

Missed that! In this case we need some kind of separate check for each since the UxGrid time varying z coordinate should only be 3 dimensions (time, node, depth).

On second thought I am also now wondering what this check is catching that the checks in the XGrid and UxGrid constructors are not, is it maybe worth just making those more robust ?

Parcels/src/parcels/_core/xgrid.py

Lines 142 to 145 in d0212c5

if "Z" in axes:

depth = ds["depth"]

if depth.ndim not in (1, 4):

raise ValueError(f"depth must be a 1D or 4D array of vertical coordinates, got {depth.ndim}D")

Parcels/src/parcels/_core/uxgrid.py

Lines 51 to 52 in d0212c5

if z.ndim not in (1, 3):

raise ValueError(f"z must be a 1D or 3D array of vertical coordinates, got {z.ndim}D")

Ah, good call! But then I'm surprised that the unit tests pass. Every 3D depth field should be throwing an Error.

@wyatt-fluidnumerics: could you dive in to see why these were not captured by this assert function? You have better feel for the code-flow on the UxGrid-side

After some digging I am pretty sure none of the unit tests caught this because they all use sigma_coordinate_lattice_dataset, which doesn't assign any dimension attributes so ds.cf.axes.get("Z") always comes back none.

import numpy as np from parcels._datasets.unstructured.generated import sigma_coordinate_lattice_dataset nx = 5 eta = np.zeros((3, nx, nx)) ds = sigma_coordinate_lattice_dataset(nx, (0.0, 4e3), (0.0, 4e3), 4, 100.0, eta=eta) print(f"zf dims: {ds['zf'].dims}, zf attrs: {ds['zf'].attrs}") print(f'ds.cf.axes.get("Z"): {ds.cf.axes.get("Z")}')

zf dims: ('time', 'zf', 'n_node'), zf attrs: {} ds.cf.axes.get("Z"): None

It is also worth noting that the new test in this PR passes with all of the assert_1D_or4D_vertical_dimension calls commented out thanks to the checks in the XGrid constructor. So it is possible that this is unneeded. The constructor checks themself are new as of #2933 and may not have been merged into the version of parcels used in #2954, which was why the original bug happened.

Hmm, even if I add

--- a/src/parcels/_datasets/unstructured/generated.py +++ b/src/parcels/_datasets/unstructured/generated.py @@ -403,4 +403,6 @@ def sigma_coordinate_lattice_dataset(nx, x_range, y_range, nz, bottom_depth, eta ds = ux.UxDataset(xr.Dataset(velocities, coords={"time": time, **vertical_coords}), uxgrid=uxgrid) ds.uxgrid.node_lon.attrs["units"] = "m" ds.uxgrid.node_lat.attrs["units"] = "m" + ds.zf.attrs["axis"] = "Z" + ds.zc.attrs["axis"] = "Z" return ds

I still get that ds.cf.axes.get("Z") is None. That seems strange?

So it is possible that this is unneeded

But the tests in xgrid are still relevant and caching 2D and 3D depth arrays, so would they not be relevant for there? Or are you only suggest we remove this in from_ugrid_conventions()?

That is strange, let me keep looking ... The tests are definitely relevant, I am just saying it is possible that all cases covered by assert_ 1D_or4D_vertical_dimensions are already covered by the constructor checks

I believe this has to do with the way UxArray internals wrap on top of xarray, the short explanation is that ds.zf.attrs['axis'] = 'Z' doesn't actually assign the attribute to the underlying ds. Instead, since UxArray data is stored in underlying xarrays, when you do .zf you are getting a new UxArray.DataArray object that wraps the underlying xarray object, so changing the zf attributes doesn't actually adjust the UxArray.DataArray that lies within ds, it just changes the wrapped copy returned by ds.zf.

This can be most easily demonstrated by running ds.zf is ds.zf on the dataset in the generator, which returns False.

For your test case the workaround is to change the attribute using ds.variables['zf'] which returns the object in ds itself rather than a wrapped copy. When you assign the z axis as an attribute that way (i.e. ds.variables'[zf'].attrs['axis'] = "Z") you get:

zf dims: ('time', 'zf', 'n_node'), zf attrs: {'axis': 'Z'} ds.cf.axes.get("Z"): ['zc', 'zf']

when running the same script from above.

Thanks, this helped! I now patched the test and the assert function in a423c45. Could you review again?

And fixing unit test to properly set Z axis

wyatt-fluidnumerics

Only a few other comments!

wyatt-fluidnumerics · 2026-10-26T20:12:04Z

    assert fset.context["c2"] == 2.0


+def test_fieldset_error_on_unsupported_depth_dimensionality():


Would it be useful to parametrize this to reach the UxGrid path too?

Added in 2531299

wyatt-fluidnumerics · 2026-10-26T21:27:52Z

+    # Assert that vertical dimensions have allowed dimensionality.
+    if (vertical_dims := ds.cf.axes.get("Z")) is not None:
+        for vertical_dim in vertical_dims:
+            ndim = np.atleast_1d(ds[vertical_dim]).ndim  # cast 0D to 1D


If the vertical dimension is time varying and loaded lazily could this not trigger an inadvertent large dask compute that instantly gets tossed?

Fixed in d34c367. From the LLM: np.atleast_1d(DataArray) converts to a numpy array, which computes lazy data. .ndim is metadata on the DataArray and costs nothing

From the LLM: np.atleast_1d(DataArray) converts to a numpy array, which computes lazy data. .ndim is metadata on the DataArray and costs nothing

Moving check on number of vertical dimensions from index_search to fi…

a02cfe0

…eldset creation

github-project-automation Bot added this to Parcels development Oct 26, 2026

github-project-automation Bot moved this to Backlog in Parcels development Oct 26, 2026

erikvansebille requested review from VeckoTheGecko and wyatt-fluidnumerics October 26, 2026 09:50

VeckoTheGecko approved these changes Oct 26, 2026

View reviewed changes

VeckoTheGecko reviewed Oct 26, 2026

View reviewed changes

Comment thread tests/test_index_search.py Outdated

Moving and renaming unit test

a423c45

wyatt-fluidnumerics approved these changes Oct 26, 2026

View reviewed changes

VeckoTheGecko reviewed Oct 26, 2026

View reviewed changes

erikvansebille added 3 commits October 26, 2026 20:16

Generalising vertical dimensionality check for uxgrid

045ae46

And fixing unit test to properly set Z axis

Merge branch 'main' into handling_vertical_dimensions_errors

966ab14

Fixing typing error

8fb6095

wyatt-fluidnumerics reviewed Oct 26, 2026

View reviewed changes

erikvansebille added 2 commits October 27, 2026 07:27

Add uxgrid test for unsupported depth dimensionality

2531299

Avoiding dask computation

d34c367

From the LLM: np.atleast_1d(DataArray) converts to a numpy array, which computes lazy data. .ndim is metadata on the DataArray and costs nothing

erikvansebille requested a review from wyatt-fluidnumerics October 27, 2026 08:37

	if "Z" in axes:
	depth = ds["depth"]
	if depth.ndim not in (1, 4):
	raise ValueError(f"depth must be a 1D or 4D array of vertical coordinates, got {depth.ndim}D")

	if z.ndim not in (1, 3):
	raise ValueError(f"z must be a 1D or 3D array of vertical coordinates, got {z.ndim}D")

		assert fset.context["c2"] == 2.0


		def test_fieldset_error_on_unsupported_depth_dimensionality():

Conversation

erikvansebille commented Oct 26, 2026

Description

Checklist

Uh oh!

VeckoTheGecko left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

wyatt-fluidnumerics left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wyatt-fluidnumerics Oct 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wyatt-fluidnumerics left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

wyatt-fluidnumerics Oct 26, 2026 •

edited

Loading