• [x] Closes #8230
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

Just to show why restoring dask arrays rather than editing encoding does not work in #8253

What happened?

There is some inconsistencies in the error raised by weighted reductions when the dimensions over which to apply the reduction are not present in all variables.

What did you expect to happen?

I think all reduction methods should have the same behaviour. I'm not sure what's the best behaviour, although I probably prefer the mean behaviour (do not raise an error when there are variables with missing dimensions).

Minimal Complete Verifiable Example

```Python import xarray as xr

ds = xr.Dataset({"foo": xr.DataArray([[1]] * 2), "bar": 1}) weighted = ds.weighted(ds["dim_0"])

weighted.mean(ds.dims) # OK weighted.std(ds.dims) # ValueError ```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

```Python

ValueError Traceback (most recent call last) Cell In[1], line 7 4 weighted = ds.weighted(ds["dim_0"]) 6 weighted.mean(ds.dims) # OK ----> 7 weighted.std(ds.dims) # ValueError

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:503, in Weighted.std(self, dim, skipna, keep_attrs) 497 def std( 498 self, 499 dim: Dims = None, 500 skipna: bool | None = None, 501 keep_attrs: bool | None = None, 502 ) -> T_Xarray: --> 503 return self._implementation( 504 self._weighted_std, dim=dim, skipna=skipna, keep_attrs=keep_attrs 505 )

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:540, in DatasetWeighted._implementation(self, func, dim, kwargs) 537 def _implementation(self, func, dim, kwargs) -> Dataset: 538 self._check_dim(dim) --> 540 return self.obj.map(func, dim=dim, **kwargs)

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/dataset.py:5964, in Dataset.map(self, func, keep_attrs, args, kwargs) 5962 if keep_attrs is None: 5963 keep_attrs = _get_keep_attrs(default=False) -> 5964 variables = { 5965 k: maybe_wrap_array(v, func(v, *args, kwargs)) 5966 for k, v in self.data_vars.items() 5967 } 5968 if keep_attrs: 5969 for k, v in variables.items():

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/dataset.py:5965, in <dictcomp>(.0) 5962 if keep_attrs is None: 5963 keep_attrs = _get_keep_attrs(default=False) 5964 variables = { -> 5965 k: maybe_wrap_array(v, func(v, args, *kwargs)) 5966 for k, v in self.data_vars.items() 5967 } 5968 if keep_attrs: 5969 for k, v in variables.items():

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:309, in Weighted._weighted_std(self, da, dim, skipna) 301 def _weighted_std( 302 self, 303 da: DataArray, 304 dim: Dims = None, 305 skipna: bool | None = None, 306 ) -> DataArray: 307 """Reduce a DataArray by a weighted std along some dimension(s).""" --> 309 return cast("DataArray", np.sqrt(self._weighted_var(da, dim, skipna)))

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:295, in Weighted._weighted_var(self, da, dim, skipna) 287 def _weighted_var( 288 self, 289 da: DataArray, 290 dim: Dims = None, 291 skipna: bool | None = None, 292 ) -> DataArray: 293 """Reduce a DataArray by a weighted var along some dimension(s).""" --> 295 sum_of_squares = self._sum_of_squares(da, dim=dim, skipna=skipna) 297 sum_of_weights = self._sum_of_weights(da, dim=dim) 299 return sum_of_squares / sum_of_weights

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:259, in Weighted._sum_of_squares(self, da, dim, skipna) 251 def _sum_of_squares( 252 self, 253 da: DataArray, 254 dim: Dims = None, 255 skipna: bool | None = None, 256 ) -> DataArray: 257 """Reduce a DataArray by a weighted sum_of_squares along some dimension(s).""" --> 259 demeaned = da - da.weighted(self.weights).mean(dim=dim) 261 return self._reduce((demeaned**2), self.weights, dim=dim, skipna=skipna)

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:483, in Weighted.mean(self, dim, skipna, keep_attrs) 477 def mean( 478 self, 479 dim: Dims = None, 480 skipna: bool | None = None, 481 keep_attrs: bool | None = None, 482 ) -> T_Xarray: --> 483 return self._implementation( 484 self._weighted_mean, dim=dim, skipna=skipna, keep_attrs=keep_attrs 485 )

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:529, in DataArrayWeighted._implementation(self, func, dim, kwargs) 528 def _implementation(self, func, dim, kwargs) -> DataArray: --> 529 self._check_dim(dim) 531 dataset = self.obj._to_temp_dataset() 532 dataset = dataset.map(func, dim=dim, **kwargs)

File ~/mambaforge/envs/xarray/lib/python3.11/site-packages/xarray/core/weighted.py:203, in Weighted._check_dim(self, dim) 201 missing_dims = set(dims) - set(self.obj.dims) - set(self.weights.dims) 202 if missing_dims: --> 203 raise ValueError( 204 f"{self.class.name} does not contain the dimensions: {missing_dims}" 205 )

ValueError: DataArrayWeighted does not contain the dimensions: {'dim_1'} ```

Anything else we need to know?

No response

Environment

What happened?

When attrs are not specified in a new Dataset/DataArray, they're changed to {} whenever they are used. Because of that, an object can slightly change silently.

What did you expect to happen?

In the example below, the tokens should be identical.

Minimal Complete Verifiable Example

```Python import xarray as xr import dask

ds = xr.Dataset({"foo": [0]}) # the assertion below would be OK if I specify attrs={} token0 = dask.base.tokenize(ds) print(ds) # this could be anything that uses attrs under the hood token1 = dask.base.tokenize(ds)

assert token0 == token1

AssertionError:

```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

No response

Anything else we need to know?

I thought we could store {} since the very beginning, adding the following here: python if attrs is None: attrs = {}

But a few tests failed, so it's probably the wrong way to fix this: ``` _________________ TestDask.test_attrs_mfdataset __________________

self = <xarray.tests.test_backends.TestDask object at 0x151cf2ef0>

def test_attrs_mfdataset(self) -> None:
    original = Dataset({"foo": ("x", np.random.randn(10))})
    with create_tmp_file() as tmp1:
        with create_tmp_file() as tmp2:
            ds1 = original.isel(x=slice(5))
            ds2 = original.isel(x=slice(5, 10))
            ds1.attrs["test1"] = "foo"
            ds2.attrs["test2"] = "bar"
            ds1.to_netcdf(tmp1)
            ds2.to_netcdf(tmp2)
            with open_mfdataset(
                [tmp1, tmp2], concat_dim="x", combine="nested"
            ) as actual:
                # presumes that attributes inherited from
                # first dataset loaded
                assert actual.test1 == ds1.test1
                # attributes from ds2 are not retained, e.g.,

              with pytest.raises(AttributeError, match=r"no attribute"):

E Failed: DID NOT RAISE <class 'AttributeError'>

/Users/mattia/MyGit/xarray/xarray/tests/test_backends.py:3576: Failed __________________ TestDask.test_open_mfdataset_attrs_file ___________________

self = <xarray.tests.test_backends.TestDask object at 0x151cf1de0>

def test_open_mfdataset_attrs_file(self) -> None:
    original = Dataset({"foo": ("x", np.random.randn(10))})
    with create_tmp_files(2) as (tmp1, tmp2):
        ds1 = original.isel(x=slice(5))
        ds2 = original.isel(x=slice(5, 10))
        ds1.attrs["test1"] = "foo"
        ds2.attrs["test2"] = "bar"
        ds1.to_netcdf(tmp1)
        ds2.to_netcdf(tmp2)
        with open_mfdataset(
            [tmp1, tmp2], concat_dim="x", combine="nested", attrs_file=tmp2
        ) as actual:
            # attributes are inherited from the master file
            assert actual.attrs["test2"] == ds2.attrs["test2"]
            # attributes from ds1 are not retained, e.g.,

          assert "test1" not in actual.attrs

E AssertionError: assert 'test1' not in {'test1': 'foo', 'test2': 'bar'} E + where {'test1': 'foo', 'test2': 'bar'} = <xarray.Dataset>\nDimensions: (x: 10)\nDimensions without coordinates: x\nData variables:\n foo (x) float64 dask.array<chunksize=(5,), meta=np.ndarray>\nAttributes:\n test1: foo\n test2: bar.attrs

/Users/mattia/MyGit/xarray/xarray/tests/test_backends.py:3594: AssertionError ________________ TestDask.test_open_mfdataset_attrs_file_path ________________

self = <xarray.tests.test_backends.TestDask object at 0x151cf1c30>

def test_open_mfdataset_attrs_file_path(self) -> None:
    original = Dataset({"foo": ("x", np.random.randn(10))})
    with create_tmp_files(2) as (tmp1, tmp2):
        tmp1 = Path(tmp1)
        tmp2 = Path(tmp2)
        ds1 = original.isel(x=slice(5))
        ds2 = original.isel(x=slice(5, 10))
        ds1.attrs["test1"] = "foo"
        ds2.attrs["test2"] = "bar"
        ds1.to_netcdf(tmp1)
        ds2.to_netcdf(tmp2)
        with open_mfdataset(
            [tmp1, tmp2], concat_dim="x", combine="nested", attrs_file=tmp2
        ) as actual:
            # attributes are inherited from the master file
            assert actual.attrs["test2"] == ds2.attrs["test2"]
            # attributes from ds1 are not retained, e.g.,

          assert "test1" not in actual.attrs

E AssertionError: assert 'test1' not in {'test1': 'foo', 'test2': 'bar'} E + where {'test1': 'foo', 'test2': 'bar'} = <xarray.Dataset>\nDimensions: (x: 10)\nDimensions without coordinates: x\nData variables:\n foo (x) float64 dask.array<chunksize=(5,), meta=np.ndarray>\nAttributes:\n test1: foo\n test2: bar.attrs

/Users/mattia/MyGit/xarray/xarray/tests/test_backends.py:3613: AssertionError ============================================================================================================================== warnings summary ============================================================================================================================== ../../miniconda3/envs/xarray/lib/python3.10/site-packages/seaborn/rcmod.py:82 /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/seaborn/rcmod.py:82: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead. if LooseVersion(mpl.version) >= "3.0":

../../miniconda3/envs/xarray/lib/python3.10/site-packages/setuptools/_distutils/version.py:346 xarray/tests/test_backends.py::TestNetCDF4Data::test_zero_dimensional_variable /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/setuptools/_distutils/version.py:346: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead. other = LooseVersion(other)

xarray/tests/test_array_api.py:10 /Users/mattia/MyGit/xarray/xarray/tests/test_array_api.py:10: UserWarning: The numpy.array_api submodule is still experimental. See NEP 47. import numpy.array_api as xp # isort:skip

xarray/tests/test_accessor_dt.py: 7 warnings xarray/tests/test_cftime_offsets.py: 5 warnings xarray/tests/test_cftimeindex.py: 64 warnings xarray/tests/test_cftimeindex_resample.py: 488 warnings xarray/tests/test_missing.py: 2 warnings /Users/mattia/MyGit/xarray/xarray/coding/times.py:365: FutureWarning: Index.ravel returning ndarray is deprecated; in a future version this will return a view on self. sample = dates.ravel()[0]

xarray/tests/test_backends.py::TestNetCDF4Data::test_zero_dimensional_variable /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/cfgrib/xarray_plugin.py:11: DeprecationWarning: distutils Version classes are deprecated. Use packaging.version instead. if LooseVersion(xr.version) <= "0.17.0":

xarray/tests/test_backends.py::TestDask::test_inline_array /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/_pytest/python.py:192: RuntimeWarning: deallocating CachingFileManager(<class 'netCDF4._netCDF4.Dataset'>, '/var/folders/_x/gdn6kyqn5d5g9j_ygdpcv5vm0000gp/T/tmp7ww27uxi/temp-2317.nc', mode='r', kwargs={'clobber': True, 'diskless': False, 'persist': False, 'format': 'NETCDF4'}), but file is not already closed. This may indicate a bug. result = testfunction(**testargs)

xarray/tests/test_backends.py::test_open_fsspec /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/fsspec/implementations/cached.py:589: ResourceWarning: unclosed file <_io.BufferedReader name='/var/folders/_x/gdn6kyqn5d5g9j_ygdpcv5vm0000gp/T/tmp8f034evp/0d56871fd8c14f69c81fd11bcd488de08bd1efce70691c6143d2a5f88be9ca84'> out[p] = open(fn, "rb").read() Enable tracemalloc to get traceback where the object was allocated. See https://docs.pytest.org/en/stable/how-to/capture-warnings.html#resource-warnings for more info.

xarray/tests/test_backends.py::test_open_fsspec /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/fsspec/implementations/cached.py:589: ResourceWarning: unclosed file <_io.BufferedReader name='/var/folders/_x/gdn6kyqn5d5g9j_ygdpcv5vm0000gp/T/tmp8f034evp/836ec38b21b701a0aae052168a2a2eab45504e6c6ba441f202e77f4f79b1a7c4'> out[p] = open(fn, "rb").read() Enable tracemalloc to get traceback where the object was allocated. See https://docs.pytest.org/en/stable/how-to/capture-warnings.html#resource-warnings for more info.

xarray/tests/test_backends.py::test_open_fsspec /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/fsspec/implementations/cached.py:589: ResourceWarning: unclosed file <_io.BufferedReader name='/var/folders/_x/gdn6kyqn5d5g9j_ygdpcv5vm0000gp/T/tmpmf5ddc2u/5b50e5f9d1df25a3c114c62d7a9dcfcd80615885572d4e2cb48894b48a393262'> out[p] = open(fn, "rb").read() Enable tracemalloc to get traceback where the object was allocated. See https://docs.pytest.org/en/stable/how-to/capture-warnings.html#resource-warnings for more info.

xarray/tests/test_backends.py::test_open_fsspec /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/fsspec/implementations/cached.py:589: ResourceWarning: unclosed file <_io.BufferedReader name='/var/folders/_x/gdn6kyqn5d5g9j_ygdpcv5vm0000gp/T/tmpmf5ddc2u/0d56871fd8c14f69c81fd11bcd488de08bd1efce70691c6143d2a5f88be9ca84'> out[p] = open(fn, "rb").read() Enable tracemalloc to get traceback where the object was allocated. See https://docs.pytest.org/en/stable/how-to/capture-warnings.html#resource-warnings for more info.

xarray/tests/test_backends.py::test_open_fsspec /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/fsspec/implementations/cached.py:589: ResourceWarning: unclosed file <_io.BufferedReader name='/var/folders/_x/gdn6kyqn5d5g9j_ygdpcv5vm0000gp/T/tmpmf5ddc2u/77aeffecc910b7c6882406131a4d36469d935a55c401298dbce90adb89d7d275'> out[p] = open(fn, "rb").read() Enable tracemalloc to get traceback where the object was allocated. See https://docs.pytest.org/en/stable/how-to/capture-warnings.html#resource-warnings for more info.

xarray/tests/test_backends.py::test_open_fsspec /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/fsspec/implementations/cached.py:589: ResourceWarning: unclosed file <_io.BufferedReader name='/var/folders/_x/gdn6kyqn5d5g9j_ygdpcv5vm0000gp/T/tmpmf5ddc2u/836ec38b21b701a0aae052168a2a2eab45504e6c6ba441f202e77f4f79b1a7c4'> out[p] = open(fn, "rb").read() Enable tracemalloc to get traceback where the object was allocated. See https://docs.pytest.org/en/stable/how-to/capture-warnings.html#resource-warnings for more info.

xarray/tests/test_calendar_ops.py: 14 warnings xarray/tests/test_cftime_offsets.py: 12 warnings xarray/tests/test_computation.py: 4 warnings /Users/mattia/MyGit/xarray/xarray/coding/cftime_offsets.py:1130: FutureWarning: Argument closed is deprecated in favor of inclusive. return pd.date_range(

xarray/tests/test_dataarray.py::TestDataArray::test_to_and_from_cdms2_sgrid xarray/tests/test_dataarray.py::TestDataArray::test_to_and_from_cdms2_sgrid xarray/tests/test_dataarray.py::TestDataArray::test_to_and_from_cdms2_sgrid xarray/tests/test_dataarray.py::TestDataArray::test_to_and_from_cdms2_sgrid /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/numpy/ma/core.py:7891: DeprecationWarning: elementwise comparison failed; this will raise an error in the future. if not np.all(xinf == filled(np.isinf(y), False)):

xarray/tests/test_dataset.py: 12 warnings xarray/tests/test_units.py: 20 warnings /Users/mattia/MyGit/xarray/xarray/core/common.py:1079: PendingDeprecationWarning: dropping variables using drop will be deprecated; using drop_vars is encouraged. cond_wdim = cond.drop(var for var in cond if dim not in cond[var].dims)

xarray/tests/test_distributed.py::test_open_mfdataset_can_open_files_with_cftime_index xarray/tests/test_distributed.py::test_open_mfdataset_can_open_files_with_cftime_index /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/tornado/ioloop.py:263: DeprecationWarning: There is no current event loop loop = asyncio.get_event_loop()

xarray/tests/test_distributed.py::test_open_mfdataset_can_open_files_with_cftime_index /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/tornado/platform/asyncio.py:326: DeprecationWarning: There is no current event loop self.old_asyncio = asyncio.get_event_loop()

xarray/tests/test_distributed.py::test_open_mfdataset_can_open_files_with_cftime_index /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/tornado/platform/asyncio.py:193: DeprecationWarning: There is no current event loop old_loop = asyncio.get_event_loop()

xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/dask/array/reductions.py:611: RuntimeWarning: All-NaN slice encountered return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)

xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] xarray/tests/test_duck_array_ops.py::test_datetime_mean[True] /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/dask/array/reductions.py:611: RuntimeWarning: All-NaN axis encountered return np.nanmin(x_chunk, axis=axis, keepdims=keepdims)

xarray/tests/test_groupby.py::test_groupby_drops_nans /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/flox/aggregate_flox.py:105: RuntimeWarning: invalid value encountered in true_divide out /= nanlen(group_idx, array, size=size, axis=axis, fill_value=0)

xarray/tests/test_plot.py::TestFacetGrid::test_facetgrid_polar xarray/tests/test_plot.py::TestFacetGrid::test_facetgrid_polar xarray/tests/test_plot.py::TestFacetGrid::test_facetgrid_polar /Users/mattia/MyGit/xarray/xarray/plot/plot.py:1478: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first. primitive = ax.pcolormesh(x, y, z, **kwargs)

xarray/tests/test_print_versions.py::test_show_versions /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/_distutils_hack/init.py:33: UserWarning: Setuptools is replacing distutils. warnings.warn("Setuptools is replacing distutils.")

xarray/tests/test_variable.py::TestVariableWithDask::test_eq_all_dtypes xarray/tests/test_variable.py::TestVariableWithDask::test_eq_all_dtypes xarray/tests/test_variable.py::TestVariableWithDask::test_eq_all_dtypes xarray/tests/test_variable.py::TestVariableWithDask::test_eq_all_dtypes /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/dask/core.py:119: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison return func(*(_execute_task(a, cache) for a in args))

xarray/tests/test_weighted.py::test_weighted_quantile_equal_weights[1-True-0.5-da2] xarray/tests/test_weighted.py::test_weighted_quantile_equal_weights[1-True-q1-da2] xarray/tests/test_weighted.py::test_weighted_quantile_equal_weights[3.14-True-0.5-da2] xarray/tests/test_weighted.py::test_weighted_quantile_equal_weights[3.14-True-q1-da2] /Users/mattia/miniconda3/envs/xarray/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1560: RuntimeWarning: All-NaN slice encountered r, k = function_base._ureduce(a,

-- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html ========================================================================================================================== short test summary info =========================================================================================================================== FAILED xarray/tests/test_backends.py::TestDask::test_attrs_mfdataset - Failed: DID NOT RAISE <class 'AttributeError'> FAILED xarray/tests/test_backends.py::TestDask::test_open_mfdataset_attrs_file - AssertionError: assert 'test1' not in {'test1': 'foo', 'test2': 'bar'} FAILED xarray/tests/test_backends.py::TestDask::test_open_mfdataset_attrs_file_path - AssertionError: assert 'test1' not in {'test1': 'foo', 'test2': 'bar'} ===================================================================================== 3 failed, 14484 passed, 1189 skipped, 211 xfailed, 65 xpassed, 671 warnings in 1223.86s (0:20:23) ====================================================================================== /Users/mattia/miniconda3/envs/xarray/lib/python3.10/multiprocessing/resource_tracker.py:224: UserWarning: resource_tracker: There appear to be 31 leaked semaphore objects to clean up at shutdown warnings.warn('resource_tracker: There appear to be %d ' ```

Environment

• [x] Closes #6970
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst ~New functions/methods are listed in api.rst~

• [x] Closes #5644
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

~New functions/methods are listed in api.rst~

Is your feature request related to a problem?

It's not related to a problem. I would find useful to have pop methods.

Describe the solution you'd like

Is it feasible to add pop methods?

For example, instead of doing this: ```python import xarray as xr

ds = xr.Dataset({"foo": None}) foo = ds["foo"] ds = ds.drop_vars("foo") I would like to be able to do this:python foo = ds.pop_vars("foo") ```

Describe alternatives you've considered

No response

Additional context

No response

When subplots are narrow (e.g., small lon range and large lat range), they overlap each other. I'm not sure, but I think that box-forced in ax.set_aspect('equal', 'box-forced') would fix the problem (#1397). However, box-forced is deprecated since Matplotlib 2.2 and I can't figure out how to fix this.

python import xarray as xr import cartopy.crs as ccrs python ds = xr.tutorial.load_dataset('air_temperature') ds = ds.sel(lon = slice(250, 300)) air = ds['air'] transform = ccrs.PlateCarree() projection = ccrs.Mercator(air.lon.values.mean(), air.lat.values.min(), air.lat.values.max()) python p = air.isel(time=[0,1]).plot(transform = transform, aspect = ds.dims['lon']/ds.dims['lat'], col = 'time', col_wrap = 1, subplot_kws = {'projection': projection}) for ax in p.axes.flat: ax.set_extent((air.lon.values.min(), air.lon.values.max(), air.lat.values.min(), air.lat.values.max()), crs = transform) ax.set_aspect('equal', 'box-forced')

Output of xr.show_versions()

• [x] Closes #7572
• [x] Tests added

~User visible changes (including notable bug fixes) are documented in whats-new.rst~ ~New functions/methods are listed in api.rst~

The latest version of netcdf-c does not allow writing a NCZarr file without the xarray's _ARRAY_DIMENSIONS attribute. Also, nczarr attributes are now lowercase.

• [x] Closes #6448
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst
• ~New functions/methods are listed in api.rst~

I think we can just use zarr to validate attributes, so we can support all types allowed by zarr.

Note that I removed the checks on the keys, as I believe we can rely on zarr for that as well. However, there is an issue with mixed types (e.g., attrs={"a": "foo", 1: "foo"}), but I think that needs to be addressed in zarr. See: https://github.com/zarr-developers/zarr-python/issues/1037

cc: @wankoelias @rabernat

What happened?

I'm getting issues using mypy with the changes to the typing of {full,zeros,ones}_like introduced in https://github.com/pydata/xarray/pull/6611.

cc: @headtr1ck

What did you expect to happen?

I think the object returned should be of the same type as the input object rather than Union[Dataset, DataArray]. For example, I don't think running mypy on the code below should return an error.

Minimal Complete Verifiable Example

```Python import xarray as xr

def test_zeros_like(da: xr.DataArray) -> xr.DataArray: return xr.zeros_like(da) ```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

Python test.py:4: error: Incompatible return value type (got "Union[Dataset, DataArray]", expected "DataArray") Found 1 error in 1 file (checked 1 source file

Anything else we need to know?

mypy 0.950 py310h6c45266_0 conda-forge mypy_extensions 0.4.3 py310h2ec42d9_5 conda-forge

Environment

What is your issue?

I noticed that the dimension order of the object returned by the latest polyval (main, unreleased) is different compared to xarray<=2022.3.0. For example, the following code returns different results. ```python import xarray as xr import numpy as np

values = np.array( [ "2021-04-01T05:25:19.000000000", "2021-04-01T05:25:29.000000000", "2021-04-01T05:25:39.000000000", "2021-04-01T05:25:49.000000000", "2021-04-01T05:25:59.000000000", "2021-04-01T05:26:09.000000000", ], dtype="datetime64[ns]", ) azimuth_time = xr.DataArray( values, name="azimuth_time", coords={"azimuth_time": values - values[0]} )

polyfit_coefficients = xr.DataArray( [ [2.33333335e-43, 1.62499999e-43, 2.79166678e-43], [-1.15316667e-30, 1.49518518e-31, 9.08833333e-31], [-2.50272583e-18, -1.23851062e-18, -2.99098229e-18], [5.83965193e-06, -1.53321770e-07, -4.84640242e-06], [4.44739216e06, 1.45053974e06, 5.29960857e06], ], dims=("degree", "axis"), coords={"axis": [0, 1, 2], "degree": [4, 3, 2, 1, 0]}, )

ds_out = xr.polyval(azimuth_time.coords["azimuth_time"], polyfit_coefficients) print(ds_out.dims) xarray v2022.3.0 ('azimuth_time', 'axis')

xarray v2022.3.1.dev103+gfc282d59 ('axis', 'azimuth_time') ``` Is this the expected behaviour? If yes, is it worth mentioning this change in what's new/breaking changes?

cc: @headtr1ck

What happened?

I'm not sure if this is a bug or an expected breaking change, but I'm not able to reproduce the results generated by polyval using a timedelta64 coordinate. The results are correct in xarray=2022.3.0, whereas they are wrong in the latest unreleased version (main, commit 6bb2b855498b5c68d7cca8cceb710365d58e604).

What did you expect to happen?

Both the stable and latest polyval functions should return the same results.

Minimal Complete Verifiable Example

```Python import xarray as xr import numpy as np

values = np.array( [ "2021-04-01T05:25:19.000000000", "2021-04-01T05:25:29.000000000", "2021-04-01T05:25:39.000000000", "2021-04-01T05:25:49.000000000", "2021-04-01T05:25:59.000000000", "2021-04-01T05:26:09.000000000", ], dtype="datetime64[ns]", ) azimuth_time = xr.DataArray( values, name="azimuth_time", coords={"azimuth_time": values - values[0]} )

polyfit_coefficients = xr.DataArray( [ [2.33333335e-43, 1.62499999e-43, 2.79166678e-43], [-1.15316667e-30, 1.49518518e-31, 9.08833333e-31], [-2.50272583e-18, -1.23851062e-18, -2.99098229e-18], [5.83965193e-06, -1.53321770e-07, -4.84640242e-06], [4.44739216e06, 1.45053974e06, 5.29960857e06], ], dims=("degree", "axis"), coords={"axis": [0, 1, 2], "degree": [4, 3, 2, 1, 0]}, )

print(xr.polyval(azimuth_time, polyfit_coefficients)) ```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

```Python

v2022.3.0 (Correct results)

<xarray.DataArray (azimuth_time: 6, axis: 3)> array([[4447392.16 , 1450539.74 , 5299608.57 ], [4505537.25588366, 1448882.82238152, 5250846.359196 ], [4563174.92026797, 1446979.12250014, 5201491.44401733], [4620298.31815291, 1444829.59596699, 5151549.377964 ], [4676900.67053846, 1442435.23739315, 5101025.78153601], [4732975.25442459, 1439797.08038974, 5049926.34223336]]) Coordinates: * azimuth_time (azimuth_time) datetime64[ns] 2021-04-01T05:25:19 ... 2021-... * axis (axis) int64 0 1 2

v2022.3.1.dev102+g6bb2b855 (Wrong results)

<xarray.DataArray (axis: 3, azimuth_time: 6)> array([[1.59620685e+30, 1.59620689e+30, 1.59620693e+30, 1.59620697e+30, 1.59620700e+30, 1.59620704e+30], [1.11164807e+30, 1.11164810e+30, 1.11164812e+30, 1.11164815e+30, 1.11164818e+30, 1.11164821e+30], [1.90975722e+30, 1.90975727e+30, 1.90975732e+30, 1.90975736e+30, 1.90975741e+30, 1.90975746e+30]]) Coordinates: * axis (axis) int64 0 1 2 * azimuth_time (azimuth_time) timedelta64[ns] 00:00:00 00:00:10 ... 00:00:50 ```

Anything else we need to know?

No response

Environment

In the example below I'm opening and concatenating two datasets using open_mfdataset. These datasets have variables with different values but overlapping coordinates. I'm concatenating along y, which is 0...4 in one dataset and 0...5 in the other. The y dimension of the resulting dataset is 0...5 which means that open_mfdataset has overwritten some values without showing any error/warning.

Is this the expected default behavior? I would expect to get at least a warning, but maybe I'm misunderstanding the default arguments.

I tried to play with the arguments, but I couldn't figure out which argument I should change to get an error in these scenarios.

MCVE Code Sample

python import xarray as xr import numpy as np

python for i in range(2): ds = xr.Dataset( {"foo": (("x", "y"), np.random.rand(4, 5 + i))}, coords={"x": np.arange(4), "y": np.arange(5 + i)}, ) print(ds) ds.to_netcdf(f"tmp{i}.nc")

<xarray.Dataset>
Dimensions:  (x: 4, y: 5)
Coordinates:
  * x        (x) int64 0 1 2 3
  * y        (y) int64 0 1 2 3 4
Data variables:
    foo      (x, y) float64 0.1271 0.6117 0.3769 0.1884 ... 0.853 0.5026 0.3762
<xarray.Dataset>
Dimensions:  (x: 4, y: 6)
Coordinates:
  * x        (x) int64 0 1 2 3
  * y        (y) int64 0 1 2 3 4 5
Data variables:
    foo      (x, y) float64 0.2841 0.6098 0.7761 0.0673 ... 0.2954 0.7212 0.3954

python DS = xr.open_mfdataset("tmp*.nc", concat_dim="y", combine="by_coords") print(DS)

<xarray.Dataset>
Dimensions:  (x: 4, y: 6)
Coordinates:
  * x        (x) int64 0 1 2 3
  * y        (y) int64 0 1 2 3 4 5
Data variables:
    foo      (x, y) float64 dask.array<chunksize=(4, 6), meta=np.ndarray>

Versions

• [x] Closes #6374
• [x] Tests added
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst
• ~New functions/methods are listed in api.rst~

What is your issue?

As part of the CZI EOSS4 grant, at B-Open we are keen on improving xarray/zarr cross-community conventions. It looks like xarray's "zarr" backend does not support Unidata NCZarr conventions while NetCDF-c>=4.8.1 supports xarray conventions.

Currently, it is possible to open a NCZarr data source using the "netcdf4" backend only. Here is an example: ```python import xarray as xr import numpy as np

ds = xr.Dataset( { "a": (("y", "x"), np.random.rand(6).reshape(2, 3)), "b": (("y", "x"), np.random.rand(6).reshape(2, 3)), }, coords={"y": [0, 1], "x": [10, 20, 30]}, ) ds.to_netcdf("file://test.nczarr#mode=nczarr")

ds_from_nczarr = xr.open_dataset("file://test.nczarr#mode=nczarr", engine="netcdf4") xr.testing.assert_identical(ds, ds_from_nczarr)

xr.open_dataset("test.nczarr", engine="zarr")

KeyError: 'Zarr object is missing the attribute _ARRAY_DIMENSIONS, which is required for xarray to determine variable dimensions.'

```

Would the community benefit from a "zarr" backend compatible with NCZarr as well? @rsignell-usgs and @rabernat suggested to discuss this within the xarray community, but I don't think a dedicated issue has been opened yet.

cc: @alexamici @aurghs @joshmoore

MCVE Code Sample

python import xarray as xr ds = xr.tutorial.load_dataset('rasm') ```python

Seasonal mean

ds_season = ds.groupby('time.season').mean() ds_season ```

<xarray.Dataset>
Dimensions:  (season: 4, x: 275, y: 205)
Coordinates:
    yc       (y, x) float64 16.53 16.78 17.02 17.27 ... 28.26 28.01 27.76 27.51
    xc       (y, x) float64 189.2 189.4 189.6 189.7 ... 17.65 17.4 17.15 16.91
  * season   (season) object 'DJF' 'JJA' 'MAM' 'SON'
Dimensions without coordinates: x, y
Data variables:
    Tair     (season, y, x) float64 nan nan nan nan ... 23.13 22.06 21.72 21.94

```python

The seasons are ordered in alphabetical order.

I want to sort them based on time.

But time was dropped, so I have to do this:

time_season = ds['time'].groupby('time.season').mean() ds_season.sortby(time_season) ```

<xarray.Dataset>
Dimensions:  (season: 4, x: 275, y: 205)
Coordinates:
    yc       (y, x) float64 16.53 16.78 17.02 17.27 ... 28.26 28.01 27.76 27.51
    xc       (y, x) float64 189.2 189.4 189.6 189.7 ... 17.65 17.4 17.15 16.91
  * season   (season) object 'SON' 'DJF' 'MAM' 'JJA'
Dimensions without coordinates: x, y
Data variables:
    Tair     (season, y, x) float64 nan nan nan nan ... 29.27 28.39 27.94 28.05

Expected Output

```python

Why does groupby drop time?

I would expect a dataset that looks like this:

ds_season['time'] = time_season ds_season ```

<xarray.Dataset>
Dimensions:  (season: 4, x: 275, y: 205)
Coordinates:
    yc       (y, x) float64 16.53 16.78 17.02 17.27 ... 28.26 28.01 27.76 27.51
    xc       (y, x) float64 189.2 189.4 189.6 189.7 ... 17.65 17.4 17.15 16.91
  * season   (season) object 'DJF' 'JJA' 'MAM' 'SON'
Dimensions without coordinates: x, y
Data variables:
    Tair     (season, y, x) float64 nan nan nan nan ... 23.13 22.06 21.72 21.94
    time     (season) object 1982-01-16 12:00:00 ... 1981-10-17 00:00:00

Problem Description

I often use groupby on time variables. When I do that, the time variable is dropped and replaced (e.g., time is replaced by season, month, year, ...). Most of the time I also want to sort the new dataset based on the original time. The example above shows why this is useful for seasons. Another example would be to sort monthly averages of a dataset that originally had daily data from Sep-2000 to Aug-2001. Why is time dropped? Does it make sense to keep it in the grouped dataset?

Output of xr.show_versions()

Code Sample, a copy-pastable example if possible

```python ds = ds.resample(time='1D',keep_attrs=True).mean()

```

Problem description

I'm downsampling in time a dataset which also contains timeless variables. I've noticed that resample adds the time dimension to the timeless variables. One workaround is: 1) Split the dataset in a timeless and a time-dependent dataset 2) Resample the time-dependent dataset 3) Merge the two datasets

This is not a big deal, but I was wondering if I'm missing some flag that avoids this behavior. If not, is it something that can be easily implemented in resample? It would be very useful for datasets with variables on staggered grids.

Output of xr.show_versions()

typing_extensions is imported in xarray/core/npcompat.py and xarray/core/utils.py.

However, typing-extensions is not a dependency in setup.cfg, and I get the following error when I install with pip install git+https://github.com/pydata/xarray using python 3.9:

python import xarray ModuleNotFoundError: No module named 'typing_extensions'

• [x] Passes pre-commit run --all-files

See: https://github.com/pydata/xarray/pull/5503#discussion_r655550526

cc: @dcherian

• [x] Closes #5495
• [x] Passes pre-commit run --all-files
• [ ] User visible changes (including notable bug fixes) are documented in whats-new.rst

• [x] Closes #3371
• [x] Tests added
• [x] Passes pre-commit run --all-files
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst
• [x] New functions/methods are listed in api.rst

• [x] Closes #5368
• [x] Tests added
• [x] Passes pre-commit run --all-files
• [ ] User visible changes (including notable bug fixes) are documented in whats-new.rst

Reduce methods such as mean drop non-numeric variables. However, there's no need to drop variables that are not actually reduced (i.e., when "reduce_dim" not in da.dims).

cc: @rcaneill

• [x] Closes #5387
• [x] Tests added
• [x] Passes pre-commit run --all-files
• [ ] User visible changes (including notable bug fixes) are documented in whats-new.rst

Lists passed to __getitem__ are converted to arrays. np.asarray infers the dtype, and for example converts a list with str and int to an array of str. Forcing dtype="O" fixes the issue.

• [x] Closes #5435
• [ ] Tests added
• [x] Passes pre-commit run --all-files
• [ ] User visible changes (including notable bug fixes) are documented in whats-new.rst

I think it does the job, although I'm not sure whether this is the best approach. A couple of questions: 1. I should probably add a test. Where is the best place? test_dask? 2. I added a chunks arguments in set_dims. I haven't updated the docstring yet... i) shape is not in the docstring either, is it because it is not meant to be public? If yes, I guess we can omit chunks as well. ii) The chunks argument I've introduced must be a dict mapping dim to chunks. Is that a problem since chunks is used elsewhere and can be slightly different (e.g., just an integer, 'auto'). iii) Should I add a check and raise an error if chunks is not a dictionary or is missing dimensions?

What happened: If I broadcast a Dataset with chunked DataArrays, the resulting DataArrays are chunked differently.

What you expected to happen: If I broadcast a dataset with 2 vectors of chunk size 1, I'd expect to get 2D arrays with chunksize (1, 1), rather than (1, N) and (M, 1).

Minimal Complete Verifiable Example:

python import xarray as xr ds = xr.Dataset({i: xr.DataArray(range(i*10), dims=f"dim_{i}").chunk(1) for i in range(1, 3)}) ds.broadcast_like(ds) <xarray.Dataset> Dimensions: (dim_1: 10, dim_2: 20) Dimensions without coordinates: dim_1, dim_2 Data variables: 1 (dim_1, dim_2) int64 dask.array<chunksize=(1, 20), meta=np.ndarray> 2 (dim_1, dim_2) int64 dask.array<chunksize=(10, 1), meta=np.ndarray>

Anything else we need to know?:

Environment:

What happened: Looks like I can't select a list of DataArrays with different name type.

What you expected to happen: If this is not a bug, consider raising a more informative error.

Minimal Complete Verifiable Example:

```python import xarray as xr from xarray import Dataset, DataArray

keys = ["foo", 1] ds = Dataset() for key in keys: ds[key] = DataArray()

ds[keys] KeyError: '1' ```

Environment:

What happened: Weighted operations don't raise an error when the dimensions passed don't exist.

What you expected to happen: This is not really a bug, but I find it a bit confusing because it's not consistent with the same "unweighted" operation.

Minimal Complete Verifiable Example:

python import xarray as xr ds = xr.tutorial.open_dataset("air_temperature") ds.weighted(xr.ones_like(ds["air"])).mean("dummy")

Environment:

• [x] Closes #5354
• [x] Tests added
• [x] Passes pre-commit run --all-files

• [x] Closes #5071
• [x] Passes pre-commit run --all-files

What happened: Looks like when I use map_blocks with functions that operate on a subset of variables, the computation is much faster if I subsample the dataset first.

What you expected to happen: In the example below, I wouldn't expect such a difference in computation time.

Minimal Complete Verifiable Example:

python import xarray as xr ds = xr.tutorial.open_dataset("air_temperature") ds["foo"] = xr.DataArray(ds["time"].values, dims="time") ds = ds.chunk(dict(time=1, lon=-1, lat=-1))

python def func(obj): return obj[["foo"]]

python %%timeit ds.map_blocks(func).compute()

17.3 s ± 179 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

```python %%timeit

Subsample the dataset before calling map_blocks

func(ds).map_blocks(func).compute() ```

5.63 s ± 175 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Anything else we need to know?:

Environment:

What happened: I think the latest pandas (1.1.0) conflicts with faceting along time dimensions.

Minimal Complete Verifiable Example:

python import xarray as xr airtemps = xr.tutorial.open_dataset("air_temperature") airtemps['air'].isel(time=slice(2)).plot(col='time') This returns the following error: KeyError: 1356998400000000000

Environment:

MCVE Code Sample

python import xarray as xr import numpy as np

```python

Create DataArray

da = xr.DataArray(np.zeros((10, 10, 4))) ```

```python

Default plot

Wrong: all of them should be 0.

da.plot(col='dim_2') ```

<xarray.plot.facetgrid.FacetGrid at 0x7f511815f3d0>

Expected Output

```python

Providing colorbar limits

Correct.

da.plot(col='dim_2', vmin=-.1, vmax=.1) ```

<xarray.plot.facetgrid.FacetGrid at 0x7f50ef14bf70>

Problem Description

If I simply use plot triggering FacetGrid to show 2D arrays which all have the same value (zero in the example above), all plots with exception of the last one are wrong. However, if I specify vmin and vmax using the same limits of the default plot, all plots are correct.

Output of xr.show_versions()

I have several repositories corresponding to different time periods of time. Each repository contains several netCDF files with different variables. Here is a simplified example: path/period1/ contains varsABC_period1.nc and varsDEF_period1.nc, and path/period2/ contains varsABC_period2.nc and varsDEF_period2.nc. All variables have the dimension T (time), but the other dimensions can be different (e.g., A has dimensions [T,X1,Y1], B has dimensions [T,X1,Y2], ...). Before version v0.11.1 a was able to easily open a dataset as follow: ```python ds=xr.open_mfdataset('path\period**.nc', concat_dim='T')

However, `open_mfdataset` now fails. It actually takes forever without returning any output, so my guess is that it's just very very slow. I've noticed that I'm able to create the dataset concatenating the files with the same variables first, then merging them:python dsABC=xr.open_mfdataset('path\period\varsABC.nc', concat_dim='T') dsDEF=xr.open_mfdataset('path\period\varsDEF.nc', concat_dim='T') ds=xr.merge([dsABC, dsDEF]) ``` I can't find what's the update in v0.11.1 that is causing this behavior. Any idea?

Output of xr.show_versions()

Code Sample

```python def cut_dataset(ds2cut, varList = ['Temp' 'S' 'Eta' 'U' 'V' 'W'], lonRange = [-180, 180], latRange = [-90, 90], depthRange = [0, float("inf")], timeRange = ['2007-09-01T00', '2008-08-31T18'], timeFreq = '1D', sampMethod = 'mean', interpC = True, saveNetCDF = False): """ Cut the dataset """

# Copy dataset
ds = ds2cut.copy(deep=True)

# Choose variables
varList_tmp = list(varList)
for varName in ds.variables:
    if all(x != 'time' for x in ds[varName].dims) or (varName=='time'):
            varList_tmp.append(varName)
toDrop = list(set(ds.variables)-set(varList_tmp))
ds = ds.drop(toDrop)

# Cut dataset
ds = ds.sel(time = slice(min(timeRange),  max(timeRange)),
            Xp1  = slice(min(lonRange),   max(lonRange)),
            Yp1  = slice(min(latRange),   max(latRange)),
            Zp1  = slice(min(depthRange), max(depthRange)))
ds = ds.sel(X    = slice(min(ds['Xp1'].values), max(ds['Xp1'].values)),
            Y    = slice(min(ds['Yp1'].values), max(ds['Yp1'].values)),
            Z    = slice(min(ds['Zp1'].values), max(ds['Zp1'].values)),
            Zu   = slice(min(ds['Zp1'].values), max(ds['Zp1'].values)),
            Zl   = slice(min(ds['Zp1'].values), max(ds['Zp1'].values)))

# Resample
if sampMethod=='snapshot':
    ds = ds.resample(time=timeFreq).first(skipna=False)
elif sampMethod=='mean':
    ds = ds.resample(time=timeFreq).mean()

# Create grid
grid = xgcm.Grid(ds,periodic=False)

# Interpolate
if interpC:
    for varName in varList:
        for dim in ds[varName].dims:
            if len(dim)>1 and dim!='time':
                ds[varName] = grid.interp(ds[varName], axis=dim[0])

# Remove useless variables
allDims = []
for varName in varList:
    for dim in ds[varName].dims:
        allDims.append(dim)
toDrop = []
for varName in ds.variables:
    if len(list(set(ds[varName].dims)-set(allDims)))>0:
        toDrop.append(varName)
ds = ds.drop(toDrop)

# Save to NetCDF
if saveNetCDF: ds.to_netcdf(saveNetCDF)

return ds, grid

3D test

ds_cut, grid_cut = cut_dataset(ds, varList = ['Eta'], latRange = [65, 65.5], depthRange = [0, 2], timeRange = ['2007-11-15T00', '2007-11-16T00'], timeFreq = '1D', sampMethod = 'mean', interpC = False, saveNetCDF = '3Dvariable.nc')

4D test

ds_cut, grid_cut = cut_dataset(ds, varList = ['Temp'], lonRange = [-30, -29.5], latRange = [65, 65.5], depthRange = [0, 2], timeRange = ['2007-11-15T00', '2007-11-16T00'], timeFreq = '1D', sampMethod = 'mean', interpC = False, saveNetCDF = '4Dvariable.nc') ```

Problem description

I'm working with a big dataset. However, most of the time I only need a small subset of data. My idea was to open and concatenate everything with open_mfdataset, and then extract subsets of data using the indexing routines. This approach works very good when I extract 3D variables (just lon, lat, and time), but it fails when I try to extract 4D variables (lon, lat, time, and depth). It doesn't actually fail, but to_netcdf takes forever. When I open a smaller dataset since the very beginning (let's say just November), then I'm also able to extract 4D variables. When I load the sub-dataset after using the indexing routines, does xarray need to read the whole original 4D variable? If yes, then I should probably change my approach and I should open subset of data since the very beginning. If no, am I doing something wrong?

Output of xr.show_versions()