What is your issue?

I was wondering if we could add some language to the docstring of xr.open_dataset regarding the chunk kwarg to make the transition for folks who have used a lot of xr.open_zarr in the past.

the current text is:

chunks (int, dict, 'auto' or None, optional) – If chunks is provided, it is used to load the new dataset into dask arrays. chunks=-1 loads the dataset with dask using a single chunk for all arrays. chunks={} loads the dataset with dask using engine preferred chunks if exposed by the backend, otherwise with a single chunk for all arrays. chunks='auto' will use dask auto chunking taking into account the engine preferred chunks. See dask chunking for more details.

I found that for opening large zarr stores, setting chunks={} reproduces the behavior of xr.open_zarr()? If this is true I think it would be great to include something like

chunks={} loads the dataset with dask using engine preferred chunks if exposed by the backend, otherwise with a single chunk for all arrays. In order to reproduce the default behavior of xr.open_zarr(...) use `xr.open_dataset(..., engine='zarr', chunks={})

to make this clear for users who have been using xr.open_zarr in the past.

Is your feature request related to a problem?

Take a simple example

python import xarray as xr da = xr.DataArray([1,2,4], dims=['x']) If I want to chunk the array I can do: python da.chunk({'x':1})

but I cant do this: python da.chunk(x=1)

TypeError: chunk() got an unexpected keyword argument 'x'

Describe the solution you'd like

I would like to be able to use the .chunk() method for dataarrays in both ways, since it is common for many xarray methods to work either way (e.g. isel()/coarsen()/etc).

@TomNicholas

Describe alternatives you've considered

No response

Additional context

No response

I found that the options for the missing_dims parameter where described inconsistently in the docstring (`"warning" vs "warn").

Not sure I found all occurrences of it here.

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

I am trying to apply an averaging function to multi year chunks of monthly model data. At the core the function performs a weighted average (and then some coordinate manipulations). I am using resample(time='1AS') and then try to map my custom function onto the data (see example below). Without actually loading the data, this step is prohibitively long in my workflow (20-30min depending on the model). Is there a way to apply this step completely lazily, like in the case where a simple non-weighted .mean() is used?

```python from dask.diagnostics import ProgressBar import xarray as xr import numpy as np

simple customized weighted mean function

def mean_func(ds): return ds.weighted(ds.weights).mean('time')

example dataset

t = xr.cftime_range(start='2000', periods=1000, freq='1AS') weights = xr.DataArray(np.random.rand(len(t)),dims=['time'], coords={'time':t}) data = xr.DataArray(np.random.rand(len(t)),dims=['time'], coords={'time':t, 'weights':weights}) ds = xr.Dataset({'data':data}).chunk({'time':1}) ds ```

Using resample with a simple mean works without any computation being triggered: python with ProgressBar(): ds.resample(time='3AS').mean('time')

But when I do the same step with my custom function, there are some computations showing up

python with ProgressBar(): ds.resample(time='3AS').map(mean_func) [########################################] | 100% Completed | 0.1s [########################################] | 100% Completed | 0.1s [########################################] | 100% Completed | 0.1s [########################################] | 100% Completed | 0.1s I am quite sure these are the same kind of computations that make my real-world workflow so slow.

I also confirmed that this not happening when I do not use resample first with ProgressBar(): mean_func(ds) this does not trigger a computation either. So this must be somehow related to resample? I would be happy to dig deeper into this, if somebody with more knowledge could point me to the right place.

Environment:

• [x] Closes #4625
• [x] Tests added
• [x] Passes isort . && black . && mypy . && flake8
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

What happened: I just updated to the newest upstream master of xarray to branch of a pull request for #4625 and noticed a strange error in my regular workflow.

I am working with a dataset ds_transformed like this:

and when I try to apply a weighted mean with python transformed_ds.o2.weighted(transformed_ds.dz_t).mean('time')

I get the following error, which I do not understand:

```

TypeError Traceback (most recent call last) <ipython-input-101-b08fde69b74b> in <module> ----> 1 test = transformed_ds.o2.weighted(transformed_ds.dz_t).mean('time')

~/code/xarray/xarray/core/common.py in weighted(self, weights) 788 """ 789 --> 790 return self._weighted_cls(self, weights) 791 792 def rolling(

~/code/xarray/xarray/core/weighted.py in init(self, obj, weights) 121 _weight_check(weights.data) 122 --> 123 self.obj = obj 124 self.weights = weights 125

TypeError: descriptor 'obj' for 'Weighted' objects doesn't apply to a 'DataArrayWeighted' object ```

Using the synthetic example from #4625 this does not show up. I am wondering if anybody has an idea what could be wrong with my dataset that would cause this error?

Anything else we need to know?:

Environment:

• [x] Closes #4541
• [x] Tests added
• [x] Passes isort . && black . && mypy . && flake8
• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

When working with large dask-array weights, this check triggers computation of the array. This affects xgcms ability to layer operations lazily

Would you be open to implement an option to skip this test, maybe with a warning displayed? Happy to submit a PR.

Code Sample, a copy-pastable example if possible

I am trying to open several ocean model data files. During the model run additional variables were written to the files. So for instance the first file will look like this:

<xarray.Dataset> Dimensions: (st_edges_ocean: 51, st_ocean: 50, time: 1, xt_ocean: 3600, yt_ocean: 2700) Coordinates: * xt_ocean (xt_ocean) float64 -279.9 -279.8 -279.7 -279.6 -279.5 ... * yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 -80.98 -80.94 ... * time (time) float64 4.401e+04 * st_ocean (st_ocean) float64 5.034 15.1 25.22 35.36 45.58 55.85 ... * st_edges_ocean (st_edges_ocean) float64 0.0 10.07 20.16 30.29 40.47 ... Data variables: jp_recycle (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> jp_reminp (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> jp_uptake (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> jo2 (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> dic_stf (time, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 2700, 3600), chunksize=(1, 2700, 3600)> o2_stf (time, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 2700, 3600), chunksize=(1, 2700, 3600)> Attributes: filename: 01210101.ocean_minibling_term_src.nc title: CM2.6_miniBling grid_type: mosaic grid_tile: 1 and the last file will look like this (with additional data variables o2_btf, dic_btf, and 'po4_btf`). ```

<xarray.Dataset> Dimensions: (st_edges_ocean: 51, st_ocean: 50, time: 1, xt_ocean: 3600, yt_ocean: 2700) Coordinates: * xt_ocean (xt_ocean) float64 -279.9 -279.8 -279.7 -279.6 -279.5 ... * yt_ocean (yt_ocean) float64 -81.11 -81.07 -81.02 -80.98 -80.94 ... * st_ocean (st_ocean) float64 5.034 15.1 25.22 35.36 45.58 55.85 ... * st_edges_ocean (st_edges_ocean) float64 0.0 10.07 20.16 30.29 40.47 ... * time (time) float64 7.25e+04 Data variables: jp_recycle (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> jp_reminp (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> jp_uptake (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> jo2 (time, st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 50, 2700, 3600), chunksize=(1, 1, 2700, 3600)> dic_stf (time, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 2700, 3600), chunksize=(1, 2700, 3600)> dic_btf (time, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 2700, 3600), chunksize=(1, 2700, 3600)> o2_stf (time, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 2700, 3600), chunksize=(1, 2700, 3600)> o2_btf (time, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 2700, 3600), chunksize=(1, 2700, 3600)> po4_btf (time, yt_ocean, xt_ocean) float64 dask.array<shape=(1, 2700, 3600), chunksize=(1, 2700, 3600)> Attributes: date: created 2014-01-08 program: time_average_netcdf.rb history: Perform time-means on all variables in 01990101.ocean_minibli... filename: 01990101.ocean_minibling_term_src.nc title: CM2.6_miniBling grid_type: mosaic grid_tile: 1 `` If I specify the additional variables to be dropped, reading all files withxarray.open_mfdataset` works like a charm. But without specifying the variables to be dropped it takes an excruciating amount of time to load.

First of all, I was wondering if there would be the possibility to display a warning if this situation occurs, suggesting to add these variables as drop_variables keyword. That would have saved me a ton of digging time.

Even better would be some way to read such datasets in a fast manner. If we could specify a fastpath option (like suggested in #1823), perhaps this could speed this task up (given that all dimensions stay the same)?

I noticed a problem the other day, when I tried to plot some very long control run data from CMIP6 on ocean.pangeo.io.

The control run of this model (CSIRO's ACCESS-ESM1-5), starts in year 101, but runs for 900 years. If I try to plot the full run as a timeseries at any point in the grid:

``` import xarray as xr import matplotlib.pyplot as plt import intake %matplotlib inline

col = intake.open_esm_datastore("https://raw.githubusercontent.com/NCAR/intake-esm-datastore/master/catalogs/pangeo-cmip6.json") cat = col.search(variable_id='o2', source_id='ACCESS-ESM1-5', experiment_id='piControl', table_id='Omon') data_dict = cat.to_dataset_dict(zarr_kwargs={'consolidated': True})

This needs a good amount of dask workers!

ds = data_dict['CMIP.CSIRO.ACCESS-ESM1-5.piControl.Omon.gn'] ds.o2.isel(i=180, j=100, lev=0).plot() ```

I get the following error:

When I restrict the timeseries to a shorter time it works fine: ds = data_dict['CMIP.CSIRO.ACCESS-ESM1-5.piControl.Omon.gn'] ds.o2.isel(i=180, j=100, lev=0, time=slice(0,500)).plot()

I assume that the internal logic tries to set the xlimit on the left side to some negative year if the full range is large. Is there a way to suppress that behavior? Or could the plotting routine default to a leftmost date of year 1 for any cftime data with gregorian calendar?

I am trying to create a function that applies a .groupby_bins operation over specified dimensions of a xarray dataset. E.g. I want to be able to sum temperture, salinity and other values grouped by oceanic oxygen concentrations. I want to be able to be flexible over which dimensions I apply the groupby_bins operation. For instance, I would like to apply it in every depth colum (resulting in an array of (x,y,time) but also over all spatial dimensions, resulting in a timeseries. I currently run into a strange error when I try the following.

Code Sample, a copy-pastable example if possible

python import xarray as xr import numpy as np import dask.array as dsa from dask.diagnostics import ProgressBar

```python def _func(data, bin_data, bins): """Group unlabeled array 'data' according to values in 'bin_data' using bins defined in 'bins' and sum all values""" labels = bins[1:] da_data = xr.DataArray(data, name='data') da_bin_data = xr.DataArray(bin_data, name='bin_data')

binned = da_data.groupby_bins(da_bin_data, bins, labels=labels,
                          include_lowest=True).sum()
return binned

def wrapper(obj, bin_obj, bins, dims): obj = obj.copy() bin_obj = bin_obj.copy() n_bins = len(bins)-1

    binned = xr.apply_ufunc(_func, obj, bin_obj, bins,
                        vectorize=True,
                        input_core_dims=[dims, dims, ['dummy']],
                        output_core_dims=[['bin_data_bins']],
                        output_dtypes=[np.float],
                        output_sizes={'bin_data_bins':n_bins},
                        dask='parallelized')
    return binned

I am showing the problem here on a sythetic example, since my current working dataset is quite big. The problem is the exact same.python

Groupby bins problem with small bins?

x_raw = np.arange(20) y_raw = np.arange(10) z_raw = np.arange(15)

x = xr.DataArray(dsa.from_array(x_raw, chunks=(-1)), dims=['x'], coords={'x':('x', x_raw)}) y = xr.DataArray(dsa.from_array(y_raw, chunks=(-1)), dims=['y'], coords={'y':('y', y_raw)}) z = xr.DataArray(dsa.from_array(z_raw, chunks=(-1)), dims=['z'], coords={'z':('z', z_raw)})

da = xr.DataArray(dsa.ones([20, 10, 15], chunks=[-1, -1, -1]), dims=['x', 'y', 'z'], coords={ 'x':x, 'y':y, 'z':z }) da ```

<xarray.DataArray 'wrapped-bb05d395159047b749ca855110244cb7' (x: 20, y: 10, z: 15)> dask.array<shape=(20, 10, 15), dtype=float64, chunksize=(20, 10, 15)> Coordinates: * x (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 * y (y) int64 0 1 2 3 4 5 6 7 8 9 * z (z) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Now I define some bins and apply the private _func on the DataArray. This works as expected. Note that the array just contains ones, hence we see 3000 in the first bin.

```python bins = np.arange(0,30,1)

apply private function on unlabled array

binned_data = _func(da.data, da.data, bins) print(binned_data.compute()) <xarray.DataArray 'data' (bin_data_bins: 29)> array([ 3000., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]) Coordinates: * bin_data_bins (bin_data_bins) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... ``` This would e.g. be an operation on a single time step. But when I now try to apply the function over the full array (core dimensions are set to all available dimensions).. I am getting a very strange error

python binned_full = wrapper(da, da, bins, dims=['x','y','z']) print(binned_full.data.compute())

```

ValueError Traceback (most recent call last) <ipython-input-13-05a5150badbf> in <module>() 1 binned_full = wrapper(da, da, bins, dims=['x','y','z']) ----> 2 print(binned_full.data.compute())

~/miniconda/envs/standard/lib/python3.6/site-packages/dask/base.py in compute(self, kwargs) 133 dask.base.compute 134 """ --> 135 (result,) = compute(self, traverse=False, kwargs) 136 return result 137

~/miniconda/envs/standard/lib/python3.6/site-packages/dask/base.py in compute(args, kwargs) 331 postcomputes = [a.dask_postcompute() if is_dask_collection(a) 332 else (None, a) for a in args] --> 333 results = get(dsk, keys, kwargs) 334 results_iter = iter(results) 335 return tuple(a if f is None else f(next(results_iter), a)

~/miniconda/envs/standard/lib/python3.6/site-packages/dask/threaded.py in get(dsk, result, cache, num_workers, kwargs) 73 results = get_async(pool.apply_async, len(pool._pool), dsk, result, 74 cache=cache, get_id=_thread_get_id, ---> 75 pack_exception=pack_exception, kwargs) 76 77 # Cleanup pools associated to dead threads

~/miniconda/envs/standard/lib/python3.6/site-packages/dask/local.py in get_async(apply_async, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, **kwargs) 519 _execute_task(task, data) # Re-execute locally 520 else: --> 521 raise_exception(exc, tb) 522 res, worker_id = loads(res_info) 523 state['cache'][key] = res

~/miniconda/envs/standard/lib/python3.6/site-packages/dask/compatibility.py in reraise(exc, tb) 58 if exc.traceback is not tb: 59 raise exc.with_traceback(tb) ---> 60 raise exc 61 62 else:

~/miniconda/envs/standard/lib/python3.6/site-packages/dask/local.py in execute_task(key, task_info, dumps, loads, get_id, pack_exception) 288 try: 289 task, data = loads(task_info) --> 290 result = _execute_task(task, data) 291 id = get_id() 292 result = dumps((result, id))

~/miniconda/envs/standard/lib/python3.6/site-packages/dask/local.py in _execute_task(arg, cache, dsk) 269 func, args = arg[0], arg[1:] 270 args2 = [_execute_task(a, cache) for a in args] --> 271 return func(*args2) 272 elif not ishashable(arg): 273 return arg

~/miniconda/envs/standard/lib/python3.6/site-packages/numpy/lib/function_base.py in call(self, args, *kwargs) 2737 vargs.extend([kwargs[_n] for _n in names]) 2738 -> 2739 return self._vectorize_call(func=func, args=vargs) 2740 2741 def _get_ufunc_and_otypes(self, func, args):

~/miniconda/envs/standard/lib/python3.6/site-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args) 2803 """Vectorized call to func over positional args.""" 2804 if self.signature is not None: -> 2805 res = self._vectorize_call_with_signature(func, args) 2806 elif not args: 2807 res = func()

~/miniconda/envs/standard/lib/python3.6/site-packages/numpy/lib/function_base.py in _vectorize_call_with_signature(self, func, args) 2844 2845 for index in np.ndindex(broadcast_shape): -> 2846 results = func((arg[index] for arg in args)) 2847 2848 n_results = len(results) if isinstance(results, tuple) else 1

<ipython-input-5-6d60d3b0b704> in _func(data, bin_data, bins) 7 8 binned = da_data.groupby_bins(da_bin_data, bins, labels=labels, ----> 9 include_lowest=True).sum() 10 return binned 11

~/Work/CODE/PYTHON/xarray/xarray/core/common.py in groupby_bins(self, group, bins, right, labels, precision, include_lowest, squeeze) 466 cut_kwargs={'right': right, 'labels': labels, 467 'precision': precision, --> 468 'include_lowest': include_lowest}) 469 470 def rolling(self, min_periods=None, center=False, **windows):

~/Work/CODE/PYTHON/xarray/xarray/core/groupby.py in init(self, obj, group, squeeze, grouper, bins, cut_kwargs) 225 226 if bins is not None: --> 227 binned = pd.cut(group.values, bins, **cut_kwargs) 228 new_dim_name = group.name + '_bins' 229 group = DataArray(binned, group.coords, name=new_dim_name)

~/miniconda/envs/standard/lib/python3.6/site-packages/pandas/core/reshape/tile.py in cut(x, bins, right, labels, retbins, precision, include_lowest) 134 precision=precision, 135 include_lowest=include_lowest, --> 136 dtype=dtype) 137 138 return _postprocess_for_cut(fac, bins, retbins, x_is_series,

~/miniconda/envs/standard/lib/python3.6/site-packages/pandas/core/reshape/tile.py in _bins_to_cuts(x, bins, right, labels, precision, include_lowest, dtype, duplicates) 227 return result, bins 228 --> 229 unique_bins = algos.unique(bins) 230 if len(unique_bins) < len(bins) and len(bins) != 2: 231 if duplicates == 'raise':

~/miniconda/envs/standard/lib/python3.6/site-packages/pandas/core/algorithms.py in unique(values) 362 363 table = htable(len(values)) --> 364 uniques = table.unique(values) 365 uniques = _reconstruct_data(uniques, dtype, original) 366

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.Int64HashTable.unique()

~/miniconda/envs/standard/lib/python3.6/site-packages/pandas/_libs/hashtable.cpython-36m-darwin.so in View.MemoryView.memoryview_cwrapper()

~/miniconda/envs/standard/lib/python3.6/site-packages/pandas/_libs/hashtable.cpython-36m-darwin.so in View.MemoryView.memoryview.cinit()

ValueError: buffer source array is read-only ``` This error only gets triggered upon computation.

Problem description

I am not sure If this is a bug or a user error on my side. I am still trying to get used to .apply_ufunc. If anybody has an idea for a workaround I would greatly appreciate it.

I am not sure if the rewrapping in xr.DataArrays in _func is actually necessary. I tried to find an equivalent functions that operates directly on dask.arrays but was not successful.

Output of xr.show_versions()

I am trying to reconstruct vertical cell depth from a z-star ocean model. This involves a few operations involving both dimensions and coordinates of a dataset like this:

<xarray.Dataset> Dimensions: (nv: 2, run: 2, st_edges_ocean: 51, st_ocean: 50, st_ocean_sub02: 10, sw_edges_ocean: 51, sw_ocean: 50, time: 240, xt_ocean: 360, xu_ocean: 360, yt_ocean: 200, yu_ocean: 200) Coordinates: area_e (yt_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> area_n (yu_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> area_t (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> area_u (yu_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> dxt (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> dxtn (yu_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> dxu (yu_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> dyt (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> dyte (yt_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> dyu (yu_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolat_c (yu_ocean, xu_ocean) float32 dask.array<shape=(200, 360), chunksize=(200, 360)> geolat_e (yt_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolat_n (yu_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolat_t (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolat_u (yu_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolon_c (yu_ocean, xu_ocean) float32 dask.array<shape=(200, 360), chunksize=(200, 360)> geolon_e (yt_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolon_n (yu_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolon_t (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> geolon_u (yu_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> ht (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> kmt (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> * nv (nv) float64 1.0 2.0 * run (run) object 'control' 'forced' * st_edges_ocean (st_edges_ocean) float64 0.0 ... 5.5e+03 * st_ocean (st_ocean) float64 5.034 ... 5.395e+03 * st_ocean_sub02 (st_ocean_sub02) float64 5.034 ... 98.62 * sw_edges_ocean (sw_edges_ocean) float64 5.034 ... 5.5e+03 * sw_ocean (sw_ocean) float64 10.07 ... 5.5e+03 * time (time) object 2181-01-16 12:00:00 ... 2200-12-16 12:00:00 tmask (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> tmask_region (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> umask (yu_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> umask_region (yu_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> wet_t (yt_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 360)> * xt_ocean (xt_ocean) float64 -279.5 ... 79.5 * xu_ocean (xu_ocean) float64 -279.0 ... 80.0 * yt_ocean (yt_ocean) float64 -81.5 -80.5 ... 89.5 * yu_ocean (yu_ocean) float64 -81.0 -80.0 ... 90.0 dst (st_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(50, 200, 360), chunksize=(50, 200, 360)> dswt (sw_ocean, yt_ocean, xt_ocean) float64 dask.array<shape=(50, 200, 360), chunksize=(50, 200, 360)> dxte (yt_ocean, xu_ocean) float64 dask.array<shape=(200, 360), chunksize=(200, 359)> dytn (yu_ocean, xt_ocean) float64 dask.array<shape=(200, 360), chunksize=(199, 360)>

The problematic step is when I assign the calculated dask.arrays to the original dataset. This happens in a function like this. python def add_vertical_spacing(ds): grid = Grid(ds) ds.coords['dst'] = calculate_ds(ds, dim='st') ds.coords['dswt'] = calculate_ds(ds, dim='sw') ds.coords['dzt'] = calculate_dz(ds['eta_t'], ds['ht'], ds['dst']) ds.coords['dzu'] = grid.min(grid.min(ds['dzt'], 'X'),'Y') return ds

This takes very long compared to a version where I assign the values as data variables: python def add_vertical_spacing(ds): grid = Grid(ds) ds.coords['dst'] = calculate_ds(ds, dim='st') ds.coords['dswt'] = calculate_ds(ds, dim='sw') ds['dzt'] = calculate_dz(ds['eta_t'], ds['ht'], ds['dst']) ds['dzu'] = grid.min(grid.min(ds['dzt'], 'X'),'Y') return ds I am not able to reproduce this problem in a smaller example yet and realize that my example is quite complex (e.g. has functions that are not shown). But I suspect that something triggers the computation of the array, when assigning a coordinate.

I have profiled my more complex code involving this function and it seems like there is a substantial increase in calls to {method 'acquire' of '_thread.lock' objects}.

Profile output of the first version (assigning coordinates)

For the second version (assigning data variables)

Does anyone have a feel for why this could happen or how I could refine my testing to get to the bottom of this?

Output of xr.show_versions()

I want to bring up an issue that has tripped up my workflow with large climate models many times. I am dealing with large data arrays of vertical cell thickness. These are 4d arrays (x, y, z, time) but I would define them as coordinates, not data_variables in the xarrays data model (e.g. they should not be multiplied by a value if a dataset is multiplied). These sort of coordinates might become more prevalent with newer ocean models like MOM6

Whenever I assign these arrays as coordinates operations on the arrays seem to trigger computation, whereas they don't if I set them up as data_variables. The example below shows this behavior. Is this a bug or done on purpose? Is there a workaround to keep these vertical thicknesses as coordinates?

``` import xarray as xr import numpy as np import dask.array as dsa

create dataset with with vertical thickness dz as data variable

data = xr.DataArray(dsa.random.random([30, 50, 200, 1000]), dims=['x','y', 'z', 't']) dz = xr.DataArray(dsa.random.random([30, 50, 200, 1000]), dims=['x','y', 'z', 't']) ds = xr.Dataset({'data':data, 'dz':dz})

another dataset with dz as coordinate

ds_new = xr.Dataset({'data':data}) ds_new.coords['dz'] = dz ```

%%time test = ds['data'] * ds['dz'] CPU times: user 1.94 ms, sys: 19.1 ms, total: 21 ms Wall time: 21.6 ms

%%time test = ds_new['data'] * ds_new['dz'] CPU times: user 17.4 s, sys: 1.98 s, total: 19.4 s Wall time: 12.5 s

Output of xr.show_versions()

Consider the following example. python import xarray as xr import numpy as np import dask.array as dsa from dask.diagnostics import ProgressBar ```

Groupby bins problem with small bins?

x_raw = np.arange(20) y_raw = np.arange(10) z_raw = np.arange(15)

x = xr.DataArray(dsa.from_array(x_raw, chunks=(-1)), dims=['x'], coords={'x':('x', x_raw)}) y = xr.DataArray(dsa.from_array(y_raw, chunks=(-1)), dims=['y'], coords={'y':('y', y_raw)}) z = xr.DataArray(dsa.from_array(z_raw, chunks=(-1)), dims=['z'], coords={'z':('z', z_raw)})

data = xr.DataArray(dsa.ones([20, 10, 15], chunks=[-1, -1, -1]), dims=['x', 'y', 'z'], coords={ 'x':x, 'y':y, 'z':z }) data <xarray.DataArray 'wrapped-bb05d395159047b749ca855110244cb7' (x: 20, y: 10, z: 15)> dask.array<shape=(20, 10, 15), dtype=float64, chunksize=(20, 10, 15)> Coordinates: * x (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 * y (y) int64 0 1 2 3 4 5 6 7 8 9 * z (z) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14

This dask array only contains ones. If I now try to apply groupby_bins with a specified array of bins (which are all below 1) it fails with a rather cryptic error. ```

```

This doesnt work

bins = np.array([0, 20, 40, 60 , 80, 100])*1e-6

binned = data.groupby_bins(data, bins).sum() binned

StopIteration Traceback (most recent call last) <ipython-input-7-dc9283bee4ea> in <module>() 2 bins = np.array([0, 20, 40, 60 , 80, 100])*1e-6 3 ----> 4 binned = data.groupby_bins(data, bins).sum() 5 binned

~/Work/CODE/PYTHON/xarray/xarray/core/common.py in wrapped_func(self, dim, axis, skipna, keep_attrs, kwargs) 20 keep_attrs=False, kwargs): 21 return self.reduce(func, dim, axis, keep_attrs=keep_attrs, ---> 22 skipna=skipna, allow_lazy=True, **kwargs) 23 else: 24 def wrapped_func(self, dim=None, axis=None, keep_attrs=False,

~/Work/CODE/PYTHON/xarray/xarray/core/groupby.py in reduce(self, func, dim, axis, keep_attrs, shortcut, kwargs) 572 def reduce_array(ar): 573 return ar.reduce(func, dim, axis, keep_attrs=keep_attrs, kwargs) --> 574 return self.apply(reduce_array, shortcut=shortcut) 575 576 ops.inject_reduce_methods(DataArrayGroupBy)

~/Work/CODE/PYTHON/xarray/xarray/core/groupby.py in apply(self, func, shortcut, kwargs) 516 applied = (maybe_wrap_array(arr, func(arr, kwargs)) 517 for arr in grouped) --> 518 return self._combine(applied, shortcut=shortcut) 519 520 def _combine(self, applied, shortcut=False):

~/Work/CODE/PYTHON/xarray/xarray/core/groupby.py in _combine(self, applied, shortcut) 520 def _combine(self, applied, shortcut=False): 521 """Recombine the applied objects like the original.""" --> 522 applied_example, applied = peek_at(applied) 523 coord, dim, positions = self._infer_concat_args(applied_example) 524 if shortcut:

~/Work/CODE/PYTHON/xarray/xarray/core/utils.py in peek_at(iterable) 114 """ 115 gen = iter(iterable) --> 116 peek = next(gen) 117 return peek, itertools.chain([peek], gen) 118

StopIteration: ```

If however the last bin includes the value 1 it runs as expected: ```

If I include a larger value at the end it works

bins = np.array([0, 20, 40, 60 , 80, 100, 1e7])*1e-6

binned = data.groupby_bins(data, bins).sum() binned <xarray.DataArray 'wrapped-bb05d395159047b749ca855110244cb7' (wrapped-bb05d395159047b749ca855110244cb7_bins: 6)> dask.array<shape=(6,), dtype=float64, chunksize=(5,)> Coordinates: * wrapped-bb05d395159047b749ca855110244cb7_bins (wrapped-bb05d395159047b749ca855110244cb7_bins) object (0.0, 2e-05] ... ```

Problem description

Is this expected behaviour? I would prefer it if it returned nan values for the bins that capture no values. It took me a bit to find out why my script using this was failing, and if this is expected behavior could a more helpful error message be considered?

Expected Output

Output of xr.show_versions()

MCVE Code Sample

We noticed a change in behavior in xarray that broke a test in xgcm.

Consider this code example: ```python import xarray as xr import numpy as np x = np.arange(5) data = xr.DataArray(np.random.rand(5), coords={'x':x }, dims=['x']) x_c = np.arange(5) + 0.5 data_c = xr.DataArray(np.random.rand(5), coords={'x_c':x_c }, dims=['x_c']) ds = xr.Dataset({'data':data, 'data_c':data_c}) del ds['data_c']

ds['x_c'] = ds['x_c'][:3] ds ```

Expected Output

In previous versions of xarray this resulted in <xarray.Dataset> Dimensions: (x: 5, x_c: 3) Coordinates: * x_c (x_c) float64 0.5 1.5 2.5 * x (x) int64 0 1 2 3 4 Data variables: data (x) float64 0.3828 0.6016 0.3603 0.414 0.35 (I tested right now with an older v0.11.3 but this works with v0.13 as well.

In the current master branch instead the coordinate gets padded with nans: <xarray.Dataset> Dimensions: (x: 5, x_c: 5) Coordinates: * x_c (x_c) float64 0.5 1.5 2.5 nan nan * x (x) int64 0 1 2 3 4 Data variables: data (x) float64 0.8369 0.4197 0.3058 0.7419 0.8126

Problem Description

We fixed the test in a new PR, but @dcherian encouraged me to submit this.

Output of xr.show_versions()

I often work with datasets that consist out of a lot of data_variables and coordinates. Often I am only concerned about a subset of variables, and for convenience drop all but a selected list of variables with a little snippet like this:

Code Sample, a copy-pastable example if possible

python def xr_keep(obj, varlist): """drop all data_vars exept the ones provided in `varlist` """ obj = obj.copy() drop_vars = [a for a in obj.data_vars if a not in varlist] return obj.drop(drop_vars)

I would love to have this functionality available as a DataArray/Dataset function. It could look something like da_slim = da.drop_all_but(['var1, 'var3']). Would this be of interest to people here? Then I could try to put in a PR.

Added example of xarray.Dataset.get method and alternative ds[var_list] syntax. - [ ] Closes #2884

Added clearer instructions on how to build xarray documentation for contribution guide - [x] Closes #2893

Code Sample, a copy-pastable example if possible

I was following the instructions to build the docs locally, using xarray/ci/requirements-py36 as test environment.

When running make html in the xarray/doc folder, I get this: ``` xarray: 0.12.1+9.gfad6d624.dirty, /Users/juliusbusecke/Work/CODE/PYTHON/xarray/xarray/init.py

Extension error: Could not import extension IPython.sphinxext.ipython_directive (exception: No module named 'IPython') make: *** [html] Error 2 (test_env) `` After installing ipython withconda install ipython`, the build works properly.

Should this line conda install -c conda-forge sphinx sphinx_rtd_theme sphinx-gallery numpydoc be updated to: conda install -c conda-forge sphinx sphinx_rtd_theme sphinx-gallery numpydoc ipython ?

Or am I overlooking something?

I am experiencing some trouble when using cartopy transformations in the plot module. I am concerned about plotting speed (plotting very high resolution maps of ocean model output takes forever). According to #657 I tried to use plot.imshow but I am getting unexpected results for the map projection. The longitude wrapping does not seem to work and as seen by the mismatch between the ocean mask data and the coastline in the example below, the projection seems to be inaccurate.

This seems to be related to the cartopy module (see the last plot which was done 'outside' of xarray).

My question is twofold, I guess

1) Is there is any other 'high speed' alternative to plot high resolution maps. 2) Since this error might not appear as drastically in all mapping scenarios, should plot.imshow display a warning, when invoked with a transformation argument, or even an error?

```python import xarray as xr %matplotlib inline import numpy as np import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt

ds = xr.open_dataset('ocean_mask.nc')

plt.figure() ax_i = plt.gca(projection=ccrs.Robinson()) ds.wet.plot.imshow(x='lonh',y='lath',ax=ax_i,transform=ccrs.PlateCarree()) ax_i.coastlines() plt.title('imshow')

plt.figure() ax_p = plt.gca(projection=ccrs.Robinson()) ds.wet.plot(x='lonh',y='lath',ax=ax_p,transform=ccrs.PlateCarree()) ax_p.coastlines() plt.title('standard plot')

plt.figure() ax = plt.gca(projection=ccrs.Robinson())

ax.imshow(ds.wet.data, transform=ccrs.PlateCarree(), extent=[ds.lonh.min().data, ds.lonh.max().data,ds.lath.min().data, ds.lath.max().data]) ax.coastlines() plt.title('cartopy imshow') ```

ocean_mask.nc.zip

This PR is meant to restore the internal plotting capabilities for objects with cftime.datetime dimensions. Based mostly on the discussions in #2164

• [x] Closes #2164
• [x] Tests added
• [x] Fully documented, including whats-new.rst for all changes and api.rst for new API

Code Sample, a copy-pastable example if possible

python path = '/work/Julius.Busecke/CM2.6_staged/CM2.6_A_V03_1PctTo2X/annual_averages' ds_ocean = xr.open_mfdataset(os.path.join(path,'ocean.*.ann.nc'), chunks={'time':1}, decode_times=False, engine='scipy') ds_ocean gives

```

ValueError Traceback (most recent call last) <ipython-input-7-278556ff881c> in <module>() 1 path = '/work/Julius.Busecke/CM2.6_staged/CM2.6_A_V03_1PctTo2X/annual_averages' ----> 2 ds_ocean = xr.open_mfdataset(os.path.join(path,'ocean.*.ann.nc'), chunks={'time':1}, decode_times=False, engine='scipy') 3 ds_ocean

~/code/miniconda/envs/standard/lib/python3.6/site-packages/xarray/backends/api.py in open_mfdataset(paths, chunks, concat_dim, compat, preprocess, engine, lock, kwargs) 503 lock = _default_lock(paths[0], engine) 504 datasets = [open_dataset(p, engine=engine, chunks=chunks or {}, lock=lock, --> 505 kwargs) for p in paths] 506 file_objs = [ds._file_obj for ds in datasets] 507

~/code/miniconda/envs/standard/lib/python3.6/site-packages/xarray/backends/api.py in <listcomp>(.0) 503 lock = _default_lock(paths[0], engine) 504 datasets = [open_dataset(p, engine=engine, chunks=chunks or {}, lock=lock, --> 505 **kwargs) for p in paths] 506 file_objs = [ds._file_obj for ds in datasets] 507

~/code/miniconda/envs/standard/lib/python3.6/site-packages/xarray/backends/api.py in open_dataset(filename_or_obj, group, decode_cf, mask_and_scale, decode_times, autoclose, concat_characters, decode_coords, engine, chunks, lock, cache, drop_variables) 283 elif engine == 'scipy': 284 store = backends.ScipyDataStore(filename_or_obj, --> 285 autoclose=autoclose) 286 elif engine == 'pydap': 287 store = backends.PydapDataStore(filename_or_obj)

~/code/miniconda/envs/standard/lib/python3.6/site-packages/xarray/backends/scipy_.py in init(self, filename_or_obj, mode, format, group, writer, mmap, autoclose) 133 filename=filename_or_obj, 134 mode=mode, mmap=mmap, version=version) --> 135 self.ds = opener() 136 self._autoclose = autoclose 137 self._isopen = True

~/code/miniconda/envs/standard/lib/python3.6/site-packages/xarray/backends/scipy_.py in _open_scipy_netcdf(filename, mode, mmap, version) 81 try: 82 return scipy.io.netcdf_file(filename, mode=mode, mmap=mmap, ---> 83 version=version) 84 except TypeError as e: # netcdf3 message is obscure in this case 85 errmsg = e.args[0]

~/code/miniconda/envs/standard/lib/python3.6/site-packages/scipy/io/netcdf.py in init(self, filename, mode, mmap, version, maskandscale) 264 265 if mode in 'ra': --> 266 self._read() 267 268 def setattr(self, attr, value):

~/code/miniconda/envs/standard/lib/python3.6/site-packages/scipy/io/netcdf.py in _read(self) 591 self._read_dim_array() 592 self._read_gatt_array() --> 593 self._read_var_array() 594 595 def _read_numrecs(self):

~/code/miniconda/envs/standard/lib/python3.6/site-packages/scipy/io/netcdf.py in _read_var_array(self) 696 # Build rec array. 697 if self.use_mmap: --> 698 rec_array = self._mm_buf[begin:begin+self._recs*self._recsize].view(dtype=dtypes) 699 rec_array.shape = (self._recs,) 700 else:

ValueError: new type not compatible with array. ```

xarray version: '0.9.6'

Problem description

I am trying to lazily read in a large number of high resolution ocean model output files. If I omit the engine='scipy' it works but takes forever. Is there a known reason why this would fail with the 'scipy' option?

I found #1313, and checked my conda environment: ``` $ conda list hdf

packages in environment at /home/Julius.Busecke/code/miniconda/envs/standard:

hdf4 4.2.12 0 conda-forge hdf5 1.8.18 1 conda-forge $ conda list netcdf

packages in environment at /home/Julius.Busecke/code/miniconda/envs/standard:

h5netcdf 0.4.2 py_0 conda-forge libnetcdf 4.4.1.1 6 conda-forge netcdf4 1.3.0 py36_0 conda-forge `` I can alsoimport netCDF4` and also load a single file using netCDF, so I am unsure if this is the same error as in #1313

I keep getting this error with some of the files for this particular model but not with others.

Any help would be greatly appreciated.

Problem description

xarrays plotting module is awesome because it detects automatically if the data is divergent or not and adjusts the colormap accordingly.

But what if I do not particularly like the default colormap choice for either data type (e.g. my current boss is not a fan of matplotlibs viridis)?

@shoyer mentioned on stackexchange that this is not currently possible.

I would be happy to submit a PR, if I can get some guidance on where to implement this.

What does everyone think about a syntax this shoule be implemented with?

Something like:

with xr.set_options(diverging_cmap='RdYlBl'): ... with xr.set_options(nondiverging_cmap='magma'): .... Or is there a better way to do this? I would personally prefer to be able to set this at the beginning of a notebook, so that it is applied to all following plots. Maybe something like the mpl.rcparams is a way to go here?

• [ ] Closes #2394 (remove if there is no corresponding issue, which should only be the case for minor changes)
• [ ] Tests passed (for all non-documentation changes)
• [ ] Tests added (for all bug fixes or enhancements)
• [ ] Fully documented, including whats-new.rst for all changes and api.rst for new API (remove if this change should not be visible to users, e.g., if it is an internal clean-up, or if this is part of a larger project that will be documented later)

I find myself editing the label on the colorbar of xarray plots often to include the units of the plotted variable. Would it be possible to extract the value of DataArray.attrs['units] (if present) and add it to the colorbar label?

Given an example DataArray with name='temperature' and attrs['units']='deg C', the current colorbar label says 'temperature'.

A 'plot_units' keyword could possibly be used like this for example:

python da.plot(plot_units=True) Then the colorbar label should be 'temperature [deg C]'.

Code Sample, a copy-pastable example if possible

```python import xarray as xr xr.version

'0.9.6'

ds = xr.open_dataset('../testing/Bianchi_o2.nc',chunks={'TIME':1}) ds

<xarray.Dataset> Dimensions: (DEPTH: 33, LATITUDE: 180, LONGITUDE: 360, TIME: 12, bnds: 2) Coordinates: * LONGITUDE (LONGITUDE) float64 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 ... * LATITUDE (LATITUDE) float64 -89.5 -88.5 -87.5 -86.5 -85.5 -84.5 -83.5 ... * DEPTH (DEPTH) float64 0.0 10.0 20.0 30.0 50.0 75.0 100.0 125.0 ... * TIME (TIME) float64 15.0 44.0 73.5 104.0 134.5 165.0 195.5 226.5 ... Dimensions without coordinates: bnds Data variables: DEPTH_bnds (DEPTH, bnds) float64 -5.0 5.0 5.0 15.0 15.0 25.0 25.0 40.0 ... TIME_bnds (TIME, bnds) float64 0.5 29.5 29.5 58.75 58.75 88.75 88.75 ... O2_LINEAR (TIME, DEPTH, LATITUDE, LONGITUDE) float64 nan nan nan nan ... Attributes: history: FERRET V5.70 (alpha) 29-Sep-11

This runs as expected

ds.isel(TIME=0).groupby_bins('O2_LINEAR', np.array([0,20,40,60,100])).max()

This crashes the kernel

ds.groupby_bins('O2_LINEAR', np.array([0,20,40,60,100])).max()

```

Problem description

I am working on ocean oxygen data and would like to compute the volume of the ocean contained within a range of concentration values.

I am trying to use groupby_bins but even with this modest size dataset (1 deg global resolution, 25 depth levels, 12 time steps) my kernel crashes every time without any error message.

I eventually want to perform this step on several TB of ocean model output, so this is concerning.

First of all I would like to ask if there is an easy way to diagnose the problem further. And secondly, are there recommendations how to compute the sum over groupby_bins for very large datasets (consisting out of dask arrays).

I noticed an undesired behavior in xarray when using xarray.open_dataset:

running the following in version 0.8.2-90-g2c7730d: import xarray as xr fid = 'dt_global_allsat_msla_uv_20140101_20140829.nc' ds = xr.open_dataset(fid) ds.time gives

<xarray.DataArray 'time' (time: 1)> array(['2013-12-31T19:00:00.000000000-0500'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 2014-01-01 Attributes: long_name: Time standard_name: time axis: T

Note the hour is 19, I also encountered files with 20. Since the time in the .nc file is given in 'days since', the expected output would be 00.

Indeed when running version 0.8.2 the output is: <xarray.DataArray 'time' (time: 1)> array(['2014-01-01T00:00:00.000000000'], dtype='datetime64[ns]') Coordinates: * time (time) datetime64[ns] 2014-01-01 Attributes: long_name: Time standard_name: time axis: T Any idea what could cause this?

Sample file used in the example can be found here: https://ufile.io/077da

I am not sure if this is a bug or user error by me.

I am basically trying to pass a SymLogNorm to all plots in a faceted plot, but it has no effect. It works when I just plot a single axis. Below I reproduced the effect with an example dataset.

Matplotlib version: 1.5.1 xarray version: 0.8.2

```python import xarray as xr import matplotlib as mpl

Here the passed 'norm=' keyword has the desired effect

da1 = xr.DataArray(np.random.exponential(size=[10,10]))

plt.figure() da1.plot()

plt.figure() da1.plot(norm=mpl.colors.SymLogNorm(0.1))

In the faceted plot this has no effect

da2 = xr.DataArray(np.random.exponential(size=[10,10,4]))

plt.figure() da.plot(x='dim_0',y='dim_1',col='dim_2',col_wrap=2)

plt.figure() da.plot(x='dim_0',y='dim_1',col='dim_2',col_wrap=3,norm=mpl.colors.SymLogNorm(0.1)) ```