Repository: xarray-contrib/cf-xarray

Issue:
wrongly shaped bounds from add_bounds method
Hey there,

was tinkering around with the `add_bounds` method and got this (from dump):
```bash
float yc_bounds(bounds, yc) ;
```

From the [CF-conventions](http://cfconventions.org/cf-conventions/cf-conventions.html#cell-boundaries), I found the following example for cell boundaries:
```bash
float A_bnds(eta, 2) ;
```
Here the shape is transposed.

Also `ncview` complains about the bounds with:
```
Error, compiled with max number of vertices for bounds var of 50!  But found a var with n=600
```

So I guess, the bounds should be transposed.

Discussion:
dcherian: Thanks @MuellerSeb . Thanks for taking the time to contribute!

Can you modify the existing tests to assert that bounds is the first dimension?

To fix the tests, you'll have to transpose the "expected" values in the tests too.
MuellerSeb: > I think you need to fix one more test with a transpose: ` def test_vertices_to_bounds():`

Was just wondering, why only one out of three asserts there fails. 
dcherian: @aulemahal can you take a look here please?
aulemahal: Last commits should have fixed the last issues in the helpers as well as adapted the new 2D bounds with this change!
MuellerSeb: LGTM. Thanks for your help!

Patch:
diff --git a/cf_xarray/accessor.py b/cf_xarray/accessor.py
index 661baca..8cabbb0 100644
--- a/cf_xarray/accessor.py
+++ b/cf_xarray/accessor.py
@@ -489,7 +489,7 @@ def _guess_bounds_dim(da, dim=None, out_dim="bounds"):
                     daXY.isel(Xbnds=1, Ybnds=0),
                 ],
                 out_dim,
-            )
+            ).transpose(..., "bounds")
         else:
             dim = dim[0]
     if dim not in da.dims:
@@ -507,7 +507,7 @@ def _guess_bounds_dim(da, dim=None, out_dim="bounds"):
     first = (bounds.isel({dim: 0}) - diff.isel({dim: 0})).assign_coords(
         {dim: da[dim][0]}
     )
-    result = xr.concat([first, bounds], dim=dim)
+    result = xr.concat([first, bounds], dim=dim).transpose(..., "bounds")
 
     return result
 
diff --git a/cf_xarray/datasets.py b/cf_xarray/datasets.py
index d41b9c8..c6a5bb4 100644
--- a/cf_xarray/datasets.py
+++ b/cf_xarray/datasets.py
@@ -156,8 +156,8 @@ romsds.coords["z_rho_dummy"] = (
 
 def _create_mollw_ds():
     # Dataset with random data on a grid that is some sort of Mollweide projection
-    XX, YY = np.mgrid[:11, :11] * 5 - 25
-    XX_bnds, YY_bnds = np.mgrid[:12, :12] * 5 - 27.5
+    YY, XX = np.mgrid[:11, :11] * 5 - 25
+    YY_bnds, XX_bnds = np.mgrid[:12, :12] * 5 - 27.5
 
     R = 50
     theta = np.arcsin(YY / (R * np.sqrt(2)))
@@ -179,7 +179,7 @@ def _create_mollw_ds():
             lon_vertices[1:, 1:],
             lon_vertices[1:, :-1],
         ),
-        axis=0,
+        axis=-1,
     )
     lat_bounds = np.stack(
         (
@@ -188,31 +188,31 @@ def _create_mollw_ds():
             lat_vertices[1:, 1:],
             lat_vertices[1:, :-1],
         ),
-        axis=0,
+        axis=-1,
     )
 
     mollwds = xr.Dataset(
         coords=dict(
             lon=xr.DataArray(
                 lon,
-                dims=("x", "y"),
+                dims=("y", "x"),
                 attrs={"units": "degrees_east", "bounds": "lon_bounds"},
             ),
             lat=xr.DataArray(
                 lat,
-                dims=("x", "y"),
+                dims=("y", "x"),
                 attrs={"units": "degrees_north", "bounds": "lat_bounds"},
             ),
         ),
         data_vars=dict(
             lon_bounds=xr.DataArray(
-                lon_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_east"}
+                lon_bounds, dims=("y", "x", "bounds"), attrs={"units": "degrees_east"}
             ),
             lat_bounds=xr.DataArray(
-                lat_bounds, dims=("bounds", "x", "y"), attrs={"units": "degrees_north"}
+                lat_bounds, dims=("y", "x", "bounds"), attrs={"units": "degrees_north"}
             ),
-            lon_vertices=xr.DataArray(lon_vertices, dims=("x_vertices", "y_vertices")),
-            lat_vertices=xr.DataArray(lat_vertices, dims=("x_vertices", "y_vertices")),
+            lon_vertices=xr.DataArray(lon_vertices, dims=("y_vertices", "x_vertices")),
+            lat_vertices=xr.DataArray(lat_vertices, dims=("y_vertices", "x_vertices")),
         ),
     )
 
diff --git a/cf_xarray/helpers.py b/cf_xarray/helpers.py
index e10aafd..b3a7c72 100644
--- a/cf_xarray/helpers.py
+++ b/cf_xarray/helpers.py
@@ -94,7 +94,7 @@ def _bounds_helper(values, n_core_dims, nbounds, order):
             top_left = values[..., -1:, :, 3]
             vertex_vals = np.block([[bot_left, bot_right], [top_left, top_right]])
         if order is None:  # We verify if the ccw version works.
-            calc_bnds = np.moveaxis(vertices_to_bounds(vertex_vals).values, 0, -1)
+            calc_bnds = vertices_to_bounds(vertex_vals).values
             order = (
                 "counterclockwise" if np.allclose(calc_bnds, values) else "clockwise"
             )
@@ -155,4 +155,6 @@ def vertices_to_bounds(
         raise ValueError(
             f"vertices format not understood. Got {vertices.dims} with shape {vertices.shape}."
         )
-    return xr.DataArray(bnd_vals, dims=out_dims[: vertices.ndim + 1])
+    return xr.DataArray(bnd_vals, dims=out_dims[: vertices.ndim + 1]).transpose(
+        ..., out_dims[0]
+    )


Test patch:
diff --git a/cf_xarray/tests/test_accessor.py b/cf_xarray/tests/test_accessor.py
index 2cc04e0..a1cb2a5 100644
--- a/cf_xarray/tests/test_accessor.py
+++ b/cf_xarray/tests/test_accessor.py
@@ -788,7 +788,9 @@ def test_add_bounds(dims):
         name = f"{dim}_bounds"
         assert name in added.coords
         assert added[dim].attrs["bounds"] == name
-        assert_allclose(added[name].reset_coords(drop=True), expected[dim])
+        assert_allclose(
+            added[name].reset_coords(drop=True), expected[dim].transpose(..., "bounds")
+        )
 
     _check_unchanged(original, ds)
 
@@ -824,7 +826,7 @@ def test_add_bounds_nd_variable():
     )
 
     actual = ds.cf.add_bounds("z", dim="x").z_bounds.reset_coords(drop=True)
-    xr.testing.assert_identical(expected, actual)
+    xr.testing.assert_identical(expected.transpose(..., "bounds"), actual)
 
     with pytest.raises(NotImplementedError):
         ds.drop_vars("x").cf.add_bounds("z", dim="x")
diff --git a/cf_xarray/tests/test_helpers.py b/cf_xarray/tests/test_helpers.py
index f14a06b..eb763f9 100644
--- a/cf_xarray/tests/test_helpers.py
+++ b/cf_xarray/tests/test_helpers.py
@@ -17,7 +17,7 @@ def test_bounds_to_vertices():
     lat_c = cfxr.bounds_to_vertices(ds.lat_bounds, bounds_dim="bounds")
     assert_array_equal(ds.lat.values + 1.25, lat_c.values[:-1])
 
-    # 2D case, CF- order
+    # 2D case
     lat_ccw = cfxr.bounds_to_vertices(
         mollwds.lat_bounds, bounds_dim="bounds", order="counterclockwise"
     )
@@ -35,13 +35,13 @@ def test_bounds_to_vertices():
     assert_equal(lon_no, lon_ccw)
 
     # Transposing the array changes the bounds direction
-    ds = mollwds.transpose("bounds", "y", "x", "y_vertices", "x_vertices")
-    lon_c = cfxr.bounds_to_vertices(
+    ds = mollwds.transpose("x", "y", "x_vertices", "y_vertices", "bounds")
+    lon_cw = cfxr.bounds_to_vertices(
         ds.lon_bounds, bounds_dim="bounds", order="clockwise"
     )
-    lon_c2 = cfxr.bounds_to_vertices(ds.lon_bounds, bounds_dim="bounds", order=None)
-    assert_equal(ds.lon_vertices, lon_c)
-    assert_equal(ds.lon_vertices, lon_c2)
+    lon_no2 = cfxr.bounds_to_vertices(ds.lon_bounds, bounds_dim="bounds", order=None)
+    assert_equal(ds.lon_vertices, lon_cw)
+    assert_equal(ds.lon_vertices, lon_no2)
 
     # Preserves dask-backed arrays
     if DaskArray is not None: