The cfgrib notebook in the documentation is broken. I guess it's related to this PR. See: https://docs.xarray.dev/en/stable/examples/ERA5-GRIB-example.html

Do you know where the in-place modification is happening? We could just copy there and fix this particular issue.

Not sure, but I'll take a look!

Hi! I was about to open a new issue about this, but looks like it's a known issue and there's a stale PR... Let me know if I can help to get this fixed!

Slightly different issue related to the latest release of netcdf-c.

Looks like nczarr attribute key changed from "_NCZARR_ARRAY" to "_nczarr_array" (see https://github.com/Unidata/netcdf-c/pull/2492), so we would have to take care of this if we want to keep supporting pure nczarr.

https://github.com/pydata/xarray/blob/6531b57f8c5cb7f3c564ff895c2e4b6573bb5521/xarray/backends/zarr.py#L202

See: https://github.com/Unidata/netcdf-c/issues/2647

The problem is that with netcdf-c=4.9.1 dimension names are lost when writing pure nczarr files. We only need it to make sure that xarray is able to read from pure nczarr.

The problem comes from libnetcdf. It bumped from 4.8.1 to 4.9.1. I'll look into it.

All set. Thanks everyone!

I'm getting errors with multi-indexes and flox. Is this expected and related to this issue, or should I open a separate issue?

```python import numpy as np

import xarray as xr

ds = xr.Dataset( dict(a=(("z",), np.ones(10))), coords=dict(b=(("z"), np.arange(2).repeat(5)), c=(("z"), np.arange(5).repeat(2))), ).set_index(bc=["b", "c"]) grouped = ds.groupby("bc")

with xr.set_options(use_flox=False): grouped.sum() # OK

with xr.set_options(use_flox=True): grouped.sum() # Error Traceback (most recent call last): File "/Users/mattia/MyGit/test.py", line 15, in <module> grouped.sum() File "/Users/mattia/MyGit/xarray/xarray/core/_reductions.py", line 2763, in sum return self._flox_reduce( File "/Users/mattia/MyGit/xarray/xarray/core/groupby.py", line 661, in _flox_reduce result = xarray_reduce( File "/Users/mattia/mambaforge/envs/sarsen_dev/lib/python3.10/site-packages/flox/xarray.py", line 373, in xarray_reduce actual[k] = v.expand_dims(missing_group_dims) File "/Users/mattia/MyGit/xarray/xarray/core/dataset.py", line 1427, in setitem self.update({key: value}) File "/Users/mattia/MyGit/xarray/xarray/core/dataset.py", line 4432, in update merge_result = dataset_update_method(self, other) File "/Users/mattia/MyGit/xarray/xarray/core/merge.py", line 1070, in dataset_update_method return merge_core( File "/Users/mattia/MyGit/xarray/xarray/core/merge.py", line 722, in merge_core aligned = deep_align( File "/Users/mattia/MyGit/xarray/xarray/core/alignment.py", line 824, in deep_align aligned = align( File "/Users/mattia/MyGit/xarray/xarray/core/alignment.py", line 761, in align aligner.align() File "/Users/mattia/MyGit/xarray/xarray/core/alignment.py", line 550, in align self.assert_unindexed_dim_sizes_equal() File "/Users/mattia/MyGit/xarray/xarray/core/alignment.py", line 450, in assert_unindexed_dim_sizes_equal raise ValueError( ValueError: cannot reindex or align along dimension 'bc' because of conflicting dimension sizes: {10, 6} (note: an index is found along that dimension with size=10) ```

Thanks - I think I might be misunderstanding how the new implementation works. I tried the following changes, but both of them return an error: python xr.polyval(values - values[0], polyfit_coefficients) Traceback (most recent call last): File "/Users/mattia/MyGit/test.py", line 31, in <module> xr.polyval(values - values[0], polyfit_coefficients) File "/Users/mattia/MyGit/xarray/xarray/core/computation.py", line 1908, in polyval coord = _ensure_numeric(coord) # type: ignore # https://github.com/python/mypy/issues/1533 ? File "/Users/mattia/MyGit/xarray/xarray/core/computation.py", line 1949, in _ensure_numeric return to_floatable(data) File "/Users/mattia/MyGit/xarray/xarray/core/computation.py", line 1939, in to_floatable x.data, ValueError: cannot include dtype 'm' in a buffer

python xr.polyval(azimuth_time.coords["azimuth_time"], polyfit_coefficients) Traceback (most recent call last): File "/Users/mattia/MyGit/test.py", line 31, in <module> xr.polyval(azimuth_time.coords["azimuth_time"], polyfit_coefficients) File "/Users/mattia/MyGit/xarray/xarray/core/computation.py", line 1908, in polyval coord = _ensure_numeric(coord) # type: ignore # https://github.com/python/mypy/issues/1533 ? File "/Users/mattia/MyGit/xarray/xarray/core/computation.py", line 1949, in _ensure_numeric return to_floatable(data) File "/Users/mattia/MyGit/xarray/xarray/core/computation.py", line 1938, in to_floatable data=datetime_to_numeric( File "/Users/mattia/MyGit/xarray/xarray/core/duck_array_ops.py", line 434, in datetime_to_numeric array = array - offset numpy.core._exceptions._UFuncBinaryResolutionError: ufunc 'subtract' cannot use operands with types dtype('<m8[ns]') and dtype('<M8[D]')

Could you also add brief updates to mention NCZarr support in the docstring for open_zarr and the user guide? In particular this paragraph should be updated:

Xarray can’t open just any zarr dataset, because xarray requires special metadata (attributes) describing the dataset dimensions and coordinates. At this time, xarray can only open zarr datasets that have been written by xarray. For implementation details, see Zarr Encoding Specification.

Documentation should be in good shape now.

Thanks for the review @shoyer! This should be good to go now.

The code now looks for NCZarr attributes if both of the following conditions are True: - _ARRAY_DIMENSIONS is missing - we are NOT in mode a or r+. As we don't write NCZarr attributes, this prevents from creating zarr using a mix of NCZarr and Xarray conventions.

I'm not sure what's the best approach with _NC* attributes. Currently, after reading the metadata useful to xarray, they are hidden or dropped. This is somewhat related to #6448.

Thanks! #6420 looks at .zarray["_NCZARR_ARRAY"]["dimrefs"] only if .zattrs["_ARRAY_ATTRIBUTE"] is missing.

Adding support for reading NCZarr in the "zarr" backend should be quite easy if xarray doesn't need to integrate the additional features in NCZarr (e.g., groups, fully qualified names, dtypes for attributes). It looks like the main difference is that the dimension names stored by xarray in .zattrs["_ARRAY_DIMENSIONS"] are stored by NCZarr in .zarray["_NCZARR_ARRAY"]["dimrefs"]. I drafted PR #6420 to explore what it would take to support reading NCZarr in xarray's "zarr" backend, and I don't think there are major changes/additions needed. (I'm experiencing issues with Windows in PR #6420. I think they need to be explored in netcdf4-python or netcdf-c though - I've added a comment in the PR)

I'm not sure whether it is better to (i) add direct support for NCZarr in xarray or (ii) just rely on the netcdf4 backend. After playing a bit with both backends, I have a few comments if option (ii) is chosen: * I would change the error raised when "_ARRAY_DIMENSIONS" is not present, suggesting to try the netcdf4 backend as well. Also, I think it's worth pointing out in the documentation or in the error message where to find information on how to open/write zarr data with the netcdf4 backend. I suspect right now it's not easy to find that information for python/xarray users. * I would consider starting a deprecation cycle for open_zarr, so it will be more clear that zarr data can be opened using various backends. * If "_ARRAY_DIMENSIONS" and "_NC*" attributes will coexist in the next version of NCZarr, the zarr backend will be able to open NCZarr but will treat "_NC*" attributes as regular attributes. I think the "zarr" backend would have to handle "_NC*" attributes (e.g., drop or hide), otherwise there can be issues when writing: TypeError: Invalid value for attr '_NCZARR_ATTR': {'types': {'Conventions': '<U1', 'title': '<U1', 'description': '<U1', 'platform': '<U1', 'references': '<U1', '_NCProperties': '<U1'}}. For serialization to netCDF files, its value must be of one of the following types: str, Number, ndarray, number, list, tuple

The errors on Windows appear to be related to the fill_value written by the "netcdf4" backend in ".zattrs", so probably needs to be addressed in netcdf4-python or netcdf-c. The fill_value is Nan rather than NaN.

@shoyer I added typing-extensions in the docs too, so you'd have to remove it from there as well: https://github.com/shoyer/xarray/blob/typing-extensions-optional/doc/getting-started-guide/installing.rst

Looks like there isn't an action that only installs xarray using pip, but dependencies are installed first using conda. Must be in the conda recipe of one of the packages specified in the CI environments?

Thanks @crusaderky! I think all your suggestions are now implemented.

So under the hood use dask broadcast_to(..., chunks=None) and users should run unify_chunks before and/or after broadcasting?

In the example above, if my target is chunksize (1, 1), ds.unify_chunks() works. Does it make any difference passing chunks=(1, 1) to dask broadcast_to rather than using chunks=None and then rechunk?

Just ran into this bug, see #5393. Hopefully no one was already working on it...

All tests are passing with dtype='O', although I'm not fully following why asarray is needed in the first place. I'll open a PR so we can start from there...

np.asarray(key, dtype='O') fixes the bug, not sure if there's any downside

I think np.asarray is converting everything to strings: https://github.com/pydata/xarray/blob/a6a1e48b57499f91db7e7c15593aadc7930020e8/xarray/core/dataset.py#L1488

Do you want to add a whatsnew @malmans2 ? Or on the next PR is fine too...

All set I think!

Got it! Thanks @keewis.

Here is the full error: ```

KeyError Traceback (most recent call last) <ipython-input-9-c00f9ae5bb67> in <module> ----> 1 airtemps['air'].isel(time=slice(2)).plot(col='time')

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/plot/plot.py in call(self, kwargs) 444 445 def call(self, kwargs): --> 446 return plot(self._da, **kwargs) 447 448 # we can't use functools.wraps here since that also modifies the name / qualname

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/plot/plot.py in plot(darray, row, col, col_wrap, ax, hue, rtol, subplot_kws, kwargs) 198 kwargs["ax"] = ax 199 --> 200 return plotfunc(darray, kwargs) 201 202

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/plot/plot.py in newplotfunc(darray, x, y, figsize, size, aspect, ax, row, col, col_wrap, xincrease, yincrease, add_colorbar, add_labels, vmin, vmax, cmap, center, robust, extend, levels, infer_intervals, colors, subplot_kws, cbar_ax, cbar_kwargs, xscale, yscale, xticks, yticks, xlim, ylim, norm, kwargs) 636 # Need the decorated plotting function 637 allargs["plotfunc"] = globals()[plotfunc.name] --> 638 return _easy_facetgrid(darray, kind="dataarray", allargs) 639 640 plt = import_matplotlib_pyplot()

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/plot/facetgrid.py in _easy_facetgrid(data, plotfunc, kind, x, y, row, col, col_wrap, sharex, sharey, aspect, size, subplot_kws, ax, figsize, kwargs) 642 643 if kind == "dataarray": --> 644 return g.map_dataarray(plotfunc, x, y, kwargs) 645 646 if kind == "dataset":

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/plot/facetgrid.py in map_dataarray(self, func, x, y, **kwargs) 263 # Get x, y labels for the first subplot 264 x, y = _infer_xy_labels( --> 265 darray=self.data.loc[self.name_dicts.flat[0]], 266 x=x, 267 y=y,

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/core/dataarray.py in getitem(self, key) 196 labels = indexing.expanded_indexer(key, self.data_array.ndim) 197 key = dict(zip(self.data_array.dims, labels)) --> 198 return self.data_array.sel(**key) 199 200 def setitem(self, key, value) -> None:

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/core/dataarray.py in sel(self, indexers, method, tolerance, drop, indexers_kwargs) 1152 method=method, 1153 tolerance=tolerance, -> 1154 indexers_kwargs, 1155 ) 1156 return self._from_temp_dataset(ds)

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/core/dataset.py in sel(self, indexers, method, tolerance, drop, **indexers_kwargs) 2100 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel") 2101 pos_indexers, new_indexes = remap_label_indexers( -> 2102 self, indexers=indexers, method=method, tolerance=tolerance 2103 ) 2104 result = self.isel(indexers=pos_indexers, drop=drop)

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/core/coordinates.py in remap_label_indexers(obj, indexers, method, tolerance, **indexers_kwargs) 395 396 pos_indexers, new_indexes = indexing.remap_label_indexers( --> 397 obj, v_indexers, method=method, tolerance=tolerance 398 ) 399 # attach indexer's coordinate to pos_indexers

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/core/indexing.py in remap_label_indexers(data_obj, indexers, method, tolerance) 268 coords_dtype = data_obj.coords[dim].dtype 269 label = maybe_cast_to_coords_dtype(label, coords_dtype) --> 270 idxr, new_idx = convert_label_indexer(index, label, dim, method, tolerance) 271 pos_indexers[dim] = idxr 272 if new_idx is not None:

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/xarray/core/indexing.py in convert_label_indexer(index, label, index_name, method, tolerance) 188 else: 189 indexer = index.get_loc( --> 190 label.item(), method=method, tolerance=tolerance 191 ) 192 elif label.dtype.kind == "b":

~/anaconda3/envs/ospy_tests/lib/python3.7/site-packages/pandas/core/indexes/datetimes.py in get_loc(self, key, method, tolerance) 620 else: 621 # unrecognized type --> 622 raise KeyError(key) 623 624 try:

KeyError: 1356998400000000000 ```

Yup, happy to do it.

Just one doubt. I think in cases where indexes[i][-1] == indexes[i+1][0], the concatenation should be consistent with the compat argument used for merge (not sure if you guys agree on this). I don't know the backend though, so the easiest thing I can think about is to run merge to trigger the exact same checks: python xr.merge([datasets[i].isel(dim=-1), datasets[i+1].isel(dim=0)], compat=compat)

Nevermind, it looks like if the check goes into _infer_concat_order_from_coords it won't affect combine_nested. So indexes[i][-1] <= indexes[i+1][0] should work.

If indexes[i] = [1, 5] and indexes[i+1] = [2, 3, 4], wouldn't indexes[i][-1] <= indexes[i+1][0] raise an error even if all indexes are different?

What about something like this? I think it would cover all possibilities, but maybe it is too expensive? python if not indexes[0].append(indexes[1:]).is_unique: raise ValueError

Got it, Thanks! Let me know if it is worth adding some checks. I'd be happy to work on it.

This fixes the problem: python plot_data = da.values kwargs = xr.plot.utils._determine_cmap_params(plot_data) if kwargs['vmin'] == kwargs['vmax']: kwargs['vmin'] -= .1 kwargs['vmax'] += .1 da.plot(col='dim_2', **kwargs)

Does it make sense to add the following somewhere in _determine_cmap_params? python if vmin == vmax: vmin -= .1 vmax += .1

I checked PR #2678 with the data that originated the issue and it fixes the problem!

I've created a little script to reproduce the problem. @TomNicholas it looks like datasets are opened correctly. The problem arises when open_mfdatasets calls _auto_combine. Indeed, _auto_combine was introduced in v0.11.1

```python import numpy as np import xarray as xr import os Tsize=100; T = np.arange(Tsize); Xsize=900; X = np.arange(Xsize); Ysize=800; Y = np.arange(Ysize) data = np.random.randn(Tsize, Xsize, Ysize) for i in range(2):

# Create 2 datasets with different variables
dsA = xr.Dataset({'A': xr.DataArray(data, coords={'T': T+i*Tsize}, dims=('T', 'X', 'Y'))})
dsB = xr.Dataset({'B': xr.DataArray(data, coords={'T': T+i*Tsize}, dims=('T', 'X', 'Y'))})

# Save datasets in one folder
dsA.to_netcdf('dsA'+str(i)+'.nc')
dsB.to_netcdf('dsB'+str(i)+'.nc')

# Save datasets in two folders
dirname='rep'+str(i)
os.mkdir(dirname)
dsA.to_netcdf(dirname+'/'+'dsA'+str(i)+'.nc')
dsB.to_netcdf(dirname+'/'+'dsB'+str(i)+'.nc')

```

Fast if netCDFs are stored in one folder:

python %%time ds_1folder = xr.open_mfdataset('*.nc', concat_dim='T')

CPU times: user 49.9 ms, sys: 5.06 ms, total: 55 ms
Wall time: 59.7 ms

Slow if netCDFs are stored in several folders:

python %%time ds_2folders = xr.open_mfdataset('rep*/*.nc', concat_dim='T')

CPU times: user 8.6 s, sys: 5.95 s, total: 14.6 s
Wall time: 10.3 s

Fast if files containing different variables are opened separately, then merged:

python %%time ds_A = xr.open_mfdataset('rep*/dsA*.nc', concat_dim='T') ds_B = xr.open_mfdataset('rep*/dsB*.nc', concat_dim='T') ds_merged = xr.merge([ds_A, ds_B])

CPU times: user 33.8 ms, sys: 3.7 ms, total: 37.5 ms
Wall time: 34.5 ms

In my previous comment I said that this would be useful for staggered grids, but then I realized that resample only operates on the time dimension. Anyway, here is my example:

```python import xarray as xr import pandas as pd import numpy as np

Create coordinates

time = pd.date_range('1/1/2018', periods=365, freq='D') space = pd.np.arange(10)

Create random variables

var_withtime1 = np.random.randn(len(time), len(space)) var_withtime2 = np.random.randn(len(time), len(space)) var_timeless1 = np.random.randn(len(space)) var_timeless2 = np.random.randn(len(space))

Create dataset

ds = xr.Dataset({'var_withtime1': (['time', 'space'], var_withtime1), 'var_withtime2': (['time', 'space'], var_withtime2), 'var_timeless1': (['space'], var_timeless1), 'var_timeless2': (['space'], var_timeless2)}, coords={'time': (['time',], time), 'space': (['space',], space)})

Standard resample: this add the time dimension to the timeless variables

ds_resampled = ds.resample(time='1M').mean()

My workaround: this does not add the time dimension to the timeless variables

ds_withtime = ds.drop([ var for var in ds.variables if not 'time' in ds[var].dims ]) ds_timeless = ds.drop([ var for var in ds.variables if 'time' in ds[var].dims ]) ds_workaround = xr.merge([ds_timeless, ds_withtime.resample(time='1M').mean()]) ```

Datasets: ```

ds <xarray.Dataset> Dimensions: (space: 10, time: 365) Coordinates: * time (time) datetime64[ns] 2018-01-01 2018-01-02 2018-01-03 ... * space (space) int64 0 1 2 3 4 5 6 7 8 9 Data variables: var_withtime1 (time, space) float64 -1.137 -0.5727 -1.287 0.8102 ... var_withtime2 (time, space) float64 1.406 0.8448 1.276 0.02579 0.5684 ... var_timeless1 (space) float64 0.02073 -2.117 -0.2891 1.735 -1.535 0.209 ... var_timeless2 (space) float64 0.4357 -0.3257 -0.8321 0.8409 0.1454 ...

ds_resampled <xarray.Dataset> Dimensions: (space: 10, time: 12) Coordinates: * time (time) datetime64[ns] 2018-01-31 2018-02-28 2018-03-31 ... * space (space) int64 0 1 2 3 4 5 6 7 8 9 Data variables: var_withtime1 (time, space) float64 0.08149 0.02121 -0.05635 0.1788 ... var_withtime2 (time, space) float64 0.08991 0.5728 0.05394 0.214 0.3523 ... var_timeless1 (time, space) float64 0.02073 -2.117 -0.2891 1.735 -1.535 ... var_timeless2 (time, space) float64 0.4357 -0.3257 -0.8321 0.8409 ...

ds_workaround <xarray.Dataset> Dimensions: (space: 10, time: 12) Coordinates: * space (space) int64 0 1 2 3 4 5 6 7 8 9 * time (time) datetime64[ns] 2018-01-31 2018-02-28 2018-03-31 ... Data variables: var_timeless1 (space) float64 0.4582 -0.6946 -0.3451 1.183 -1.14 0.1849 ... var_timeless2 (space) float64 1.658 -0.1719 -0.2202 -0.1789 -1.247 ... var_withtime1 (time, space) float64 -0.3901 0.3725 0.02935 -0.1315 ... var_withtime2 (time, space) float64 0.07145 -0.08536 0.07049 0.1025 ... ```

Alright, I found the problem. I'm loading several variables from different files. All the variables have 1464 snapshots. However, one of the 3D variables has just one snapshot at a different time (I found a bag in my bash script to re-organize the raw data). When I load my dataset using .open_mfdataset, the time dimension has an extra snapshot (length is 1465). However, xarray doesn't like it and when I run functions such as to_netcdf it takes forever (no error). Thanks @fujiisoup for the help!

I forgot to mention that I'm getting this warning: /home/idies/anaconda3/lib/python3.5/site-packages/dask/core.py:306: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison elif type_arg is type(key) and arg == key:

However, I don't think it is relevant since I get the same warning when I'm able to run .to_netcdf() on the 3D variable.

The problem occurs when I run the very last line, which is to_netcdf(). Right before, the dataset looks like this: python <xarray.Dataset> Dimensions: (X: 10, Y: 25, Z: 1, time: 2) Coordinates: * time (time) datetime64[ns] 2007-11-15 2007-11-16 * Z (Z) float64 1.0 * X (X) float64 -29.94 -29.89 -29.85 -29.81 -29.76 -29.72 -29.67 ... * Y (Y) float64 65.01 65.03 65.05 65.07 65.09 65.11 65.13 65.15 ... Data variables: drF (time, Z) float64 2.0 2.0 dxF (time, Y, X) float64 2.066e+03 2.066e+03 2.066e+03 2.066e+03 ... dyF (time, Y, X) float64 2.123e+03 2.123e+03 2.123e+03 2.123e+03 ... rA (time, Y, X) float64 4.386e+06 4.386e+06 4.386e+06 4.386e+06 ... fCori (time, Y, X) float64 0.0001322 0.0001322 0.0001322 0.0001322 ... R_low (time, Y, X) float64 -2.001e+03 -1.989e+03 -1.973e+03 ... Ro_surf (time, Y, X) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... Depth (time, Y, X) float64 2.001e+03 1.989e+03 1.973e+03 1.963e+03 ... HFacC (time, Z, Y, X) float64 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ... Temp (time, Z, Y, X) float64 dask.array<shape=(2, 1, 25, 10), chunksize=(1, 1, 25, 10)> This is a dask array, right?

I have the same issue if I don't copy the dataset.

Here are the coordinates of my dataset: python <xarray.Dataset> Dimensions: (X: 960, Xp1: 961, Y: 880, Yp1: 881, Z: 216, Zl: 216, Zp1: 217, Zu: 216, time: 1465) Coordinates: * Z (Z) float64 1.0 3.5 7.0 11.5 17.0 23.5 31.0 39.5 49.0 59.5 ... * Zp1 (Zp1) float64 0.0 2.0 5.0 9.0 14.0 20.0 27.0 35.0 44.0 54.0 ... * Zu (Zu) float64 2.0 5.0 9.0 14.0 20.0 27.0 35.0 44.0 54.0 65.0 ... * Zl (Zl) float64 0.0 2.0 5.0 9.0 14.0 20.0 27.0 35.0 44.0 54.0 ... * X (X) float64 -46.92 -46.83 -46.74 -46.65 -46.57 -46.48 -46.4 ... * Y (Y) float64 56.81 56.85 56.89 56.93 56.96 57.0 57.04 57.08 ... * Xp1 (Xp1) float64 -46.96 -46.87 -46.78 -46.7 -46.61 -46.53 ... * Yp1 (Yp1) float64 56.79 56.83 56.87 56.91 56.95 56.98 57.02 ... * time (time) datetime64[ns] 2007-09-01 2007-09-01T06:00:00 ... I don't think the horizontal coordinates are the problem because it works fine when I use the same function on 3D variables. I'm also attaching the function that I use to open the dataset, just in case is helpful: ```python def load_dataset(): """ Load the whole dataset """

# Import grid and fields separately, then merge
gridpath = '/home/idies/workspace/OceanCirculation/exp_ASR/grid_glued.nc'
fldspath = '/home/idies/workspace/OceanCirculation/exp_ASR/result_*/output_glued/*.*_glued.nc'
gridset = xr.open_dataset(gridpath,
                          drop_variables = ['XU','YU','XV','YV','RC','RF','RU','RL'])
fldsset = xr.open_mfdataset(fldspath,
                            concat_dim     = 'T',
                            drop_variables = ['diag_levels','iter'])
ds = xr.merge([gridset, fldsset])

# Adjust dimensions creating conflicts
ds = ds.rename({'Z': 'Ztmp'})
ds = ds.rename({'T': 'time', 'Ztmp': 'Z', 'Zmd000216': 'Z'})
ds = ds.squeeze('Zd000001')
for dim in ['Z','Zp1', 'Zu','Zl']:
    ds[dim].values   = np.fabs(ds[dim].values)
    ds[dim].attrs.update({'positive': 'down'})

# Create horizontal vectors (remove zeros due to exch2)
ds['X'].values   = ds.XC.where((ds.XC!=0) & (ds.YC!=0)).mean(dim='Y',   skipna=True)
ds['Xp1'].values = ds.XG.where((ds.XG!=0) & (ds.YG!=0)).mean(dim='Yp1', skipna=True)
ds['Y'].values   = ds.YC.where((ds.XC!=0) & (ds.YC!=0)).mean(dim='X',   skipna=True)
ds['Yp1'].values = ds.YG.where((ds.XG!=0) & (ds.YG!=0)).mean(dim='Xp1', skipna=True)
ds = ds.drop(['XC','YC','XG','YG'])

# Create xgcm grid
ds['Z'].attrs.update({'axis': 'Z'})
ds['X'].attrs.update({'axis': 'X'})
ds['Y'].attrs.update({'axis': 'Y'})
for dim in ['Zp1','Zu','Zl','Xp1','Yp1']:
    if min(ds[dim].values)<min(ds[dim[0]].values):
        ds[dim].attrs.update({'axis': dim[0], 'c_grid_axis_shift': -0.5})
    elif min(ds[dim].values)>min(ds[dim[0]].values):
        ds[dim].attrs.update({'axis': dim[0], 'c_grid_axis_shift': +0.5})
grid = xgcm.Grid(ds,periodic=False)

return ds, grid

``` I think somewhere I trigger the loading of the whole dataset. Otherwise, I don't understand why it works when I open just one month instead of the whole year.