What is your issue?

In #4179 / #4907, the xarray policy around minimum supported version of dependencies was changed, with the reasoning that the previous policy (based on NEP-29) was too aggressive. Ironically, this caused xarray to drop Python 3.8 on Jan 26th (#7461), 3 months before what NEP-29 recommends (Apr 14th). This is hard to defend - and in fact it sparked discontent (see late comments in #7461).

Regardless of what policy xarray decides to use internally, it should never be more aggressive than NEP-29. The xarray documentation is also incorrect, as it states "Python: 24 months (NEP-29)" which is not, in fact, in NEP-29.

It should be relatively straightforward to allow slicing on coordinates that are not backed by an IndexVariable, or in other words coordinates that are on a dimension with a different name, as long as they are 1-dimensional (unsure about the multidimensional case).

E.g. given this array: ``` a = xarray.DataArray( [10, 20, 30], dims=['country'], coords={ 'country': ['US', 'Germany', 'France'], 'currency': ('country', ['USD', 'EUR', 'EUR']) })

<xarray.DataArray (country: 3)> array([10, 20, 30]) Coordinates: * country (country) <U7 'US' 'Germany' 'France' currency (country) <U3 'USD' 'EUR' 'EUR' ```

This is currently not possible: ``` a.sel(currency='EUR')

ValueError: dimensions or multi-index levels ['currency'] do not exist ```

It should be interpreted as a shorthand for: ``` a.sel(country=a.currency == 'EUR')

<xarray.DataArray (country: 2)> array([20, 30]) Coordinates: * country (country) <U7 'Germany' 'France' currency (country) <U3 'EUR' 'EUR' ```

unstack() should be smart enough to recognise that all labels in a coord are strings, and convert them to numpy strings. This is particularly relevant e.g. if you want to dump the xarray to netcdf and then read it with a non-python library.

``` python import xarray

a = xarray.DataArray([[1,2],[3,4]], dims=['x', 'y'], coords={'x': ['x1', 'x2'], 'y': ['y1', 'y2']}) a ```

<xarray.DataArray (x: 2, y: 2)> array([[1, 2], [3, 4]]) Coordinates: * y (y) <U2 'y1' 'y2' * x (x) <U2 'x1' 'x2'

python a.stack(s=['x', 'y']).unstack('s')

<xarray.DataArray (x: 2, y: 2)> array([[1, 2], [3, 4]]) Coordinates: * x (x) object 'x1' 'x2' * y (y) object 'y1' 'y2'

The documentation for set_options laconically says:

arithmetic_join: DataArray/Dataset alignment in binary operations. Default: 'inner'.

leaving the user wonder what the other options are. Also, the set_options code does not make any kind of domain check on the possible values. By scanning the code I gathered that the valid values (and their meanings) should be the same as align(join=...), but I'd like confirmation on that...

sphinx-build emits a ton of errors that need to be polished out:

https://readthedocs.org/projects/xray/builds/ -> latest -> open last step

Options for the long term: - Change the "Docs" azure pipelines job to crash if there are new failures. From past experience though, this should come together with a sensible way to whitelist errors that can't be fixed. This will severely slow down development as PRs will systematically fail on such a check. - Add a task in the release process where, immediately before closing a release, the maintainer needs to manually go through the sphinx-build log and fix any new issues. This would be a major extra piece of work for the maintainer.

I am honestly not excited by either of the above. Alternative suggestions are welcome.

Iterative improvement on #3276 @dcherian

map_blocks crashes with ImportError if dask isn't installed, even if it's legal to run it on a DataArray/Dataset without any dask variables. This forces writers of extension libraries to either not use map_blocks, add dask as a strict requirement, or write a switch in their own code.

Please change the code so that it works without dask (you'll need to write a stub of dask.is_dask_collection that always returns False) and add relevant tests to be triggered in our py36-bare-minimum CI environment.

This is an offshoot of #3358.

The readthedocs CI has a bad habit of failing even after the Azure Pipelines job "Docs" has succeeded.

After major changes that impact the documentation, and before merging everything into master, it would be advisable to explicitly verify that RTD builds correctly.

So far I tried to 1. create my own readthedocs project, https://readthedocs.org/projects/crusaderky-xarray/ 2. point it to my fork https://github.com/crusaderky/xarray/ 3. enable build for the branch I want to merge

This is currently failing because of an issue with versioneer, which incorrectly sets xarray.__version__ to 0+untagged.111.g6d60700. This in turn causes a failure in a minimum version check in pandas.DataFrame.to_xarray() on pandas>=0.25.

In the master RTD project https://readthedocs.org/projects/xray/, I can instead read xarray: 0.13.0+20.gdd2b803a.

So far the only workaround I could find was to downgrade pandas to 0.24 in ci/requirements/doc.yml.

I think with the multiple recent breakages we've just had due to dependency upgrades, we should push out a patch release with some haste.

Please comment/add/object

Must have

• [x] numpy 1.18 support #3409
• [x] pseudonetcdf 3.1.0 support #3409, #3420
• [x] require cftime != 1.0.4 #3463
• [x] groupby reduce regression fix #3403
• [x] pandas master support #3440

Nice to have

• [x] ellipsis (...) work #1081, #3414, #3418, #3421, #3423, #3424
• [x] HTML repr #3425 (really mouth-watering, but I'm unsure about how far it is from completion)
• [x] groupby drop nan groups #3406
• [x] deprecate allow_lazy #3435
• [x] __dask_tokenize__ #3446
• [x] dask name equality #3453
• [x] Leave empty slot when not using accessors #3531

```

dict(xarray.DataArray([1, 2]).chunk().dask_graph()) { ('xarray-<this-array>-7e885b8e329090da3fe58d4483c0cf8b', 0): (<function dask.array.core.getter(a, b, asarray=True, lock=None)>, 'xarray-<this-array>-7e885b8e329090da3fe58d4483c0cf8b', (slice(0, 2, None),)), 'xarray-<this-array>-7e885b8e329090da3fe58d4483c0cf8b': ImplicitToExplicitIndexingAdapter(array=NumpyIndexingAdapter(array=array([1, 2]))) } There is no reason why this should be any more complicated than da.from_array: dict(da.from_array(np.array([1, 2]), chunks=2).dask_graph()) { ('array-de932becc43e72c010bc91ffefe42af1', 0): (<function dask.array.core.getter(a, b, asarray=True, lock=None)>, 'array-original-de932becc43e72c010bc91ffefe42af1', (slice(0, 2, None),)), 'array-original-de932becc43e72c010bc91ffefe42af1': array([1, 2]) } ``` da.from_array itself should be simplified - see twin issue https://github.com/dask/dask/issues/3556

@dcherian The new page http://xarray.pydata.org/en/stable/howdoi.html (#3357) is somewhat painful to read on readthedocs. The table goes out of the screen and one is forced to scroll left and right non stop.

Maybe a better alternative could be with Sphinx definitions syntax (which allows for automatic reflowing)?

rst How do I ... ============ Add variables from other datasets to my dataset? :py:meth:`Dataset.merge` (that's a 4 spaces indent)

In CI there are:

• tests for all the latest versions of all libraries, mandatory and optional (py36, py37, py37-windows)
• tests for the minimum versions of the mandatory libraries only (py35-min)

There are no tests for legacy versions of the optional libraries.

Today I tried downgrading dask in the py37 environment to dask=1.1.2, which is 6 months old...

...it's a bloodbath. 383 errors of the most diverse kind.

In the codebase I found mentions to much older minimum versions: installing.rst mentions dask >=0.16.1, and Dataset.chunk() even asks for dask>=0.9.

It think we should add CI tests for old versions of the optional dependencies. What policy should we adopt when we find an incompatibility? How old a library should be not to bother fixing bugs and just require a newer version? I personally would go for an aggressive 6 months worth' of backwards compatibility; less if the time it takes to fix the issues is excessive. The tests should run on py36 because py35 builds are becoming very scarce in anaconda.

This has the outlook of being an exercise in extreme frustration. I'm afraid I personally hold zero interest towards packages older than the latest available in the anaconda official repo, so I'm not volunteering for this one (sorry).

I'd like to hear other people's opinions and/or offers of self-immolation... :)

On Saturday night, all Python 3.6 CI builds started failing. Python 3.7 is unaffected. See https://dev.azure.com/xarray/xarray/_build/results?buildId=362&view=logs

MacOSX py36: UnsatisfiableError: The following specifications were found to be in conflict: - pynio - python=3.6 - rasterio

Linux py36: UnsatisfiableError: The following specifications were found to be in conflict: - cfgrib[version='>=0.9.2'] - h5netcdf - pynio

Azure pipelines is not being triggered by PRs this morning. See https://github.com/pydata/xarray/pull/3358 and https://github.com/pydata/xarray/pull/3365.

Last run was 12 hours ago.

3117 by @nvictus introduces support for sparse in plain xarray.

dask already supports it.

Running with: - xarray git head - dask 2.2.0 - numpy 1.16.4 - sparse 0.7.0 - NUMPY_EXPERIMENTAL_ARRAY_FUNCTION=1

```python

import numpy, sparse, xarray, dask.array s = sparse.COO(numpy.array([1, 2]))
da1 = dask.array.from_array(s) da1._meta <COO: shape=(0,), dtype=int64, nnz=0, fill_value=0> da1.compute() <COO: shape=(2,), dtype=int64, nnz=2, fill_value=0> da2 = xarray.DataArray(s).chunk().data da2._meta
array([], dtype=int64) # Wrong da2.compute() RuntimeError: Cannot convert a sparse array to dense automatically. To manually densify, use the todense method. ```

In my project I'm struggling with weighted sums of 2000-4000 dask-based xarrays. The time to reach the final dask-based array, the size of the final dask dict, and the time to compute the actual result are horrendous.

So I wrote the below which - as laborious as it may look - gives a performance boost nothing short of miraculous. At the bottom you'll find some benchmarks as well.

https://gist.github.com/crusaderky/62832a5ffc72ccb3e0954021b0996fdf

In my project, this deflated the size of the final dask dict from 5.2 million keys to 3.3 million and cut a 30% from the time required to define it.

I think it's generic enough to be a good addition to the core xarray module. Impressions?

I, like many others, have irreversibly fallen in love with black. Can we apply it to the existing codebase and as an enforced CI test? The only (big) problem is that developers will need to manually apply it to any open branches and then merge from master - and even then, merging likely won't be trivial. How did the dask project tackle the issue?

The 'Docs' ci project is failing because http://naciscdn.org is unresponsive:

https://dev.azure.com/xarray/xarray/_build/results?buildId=408&view=logs&jobId=7e620c85-24a8-5ffa-8b1f-642bc9b1fc36

Excerpt: ``` /usr/share/miniconda/envs/xarray-tests/lib/python3.7/site-packages/cartopy/io/init.py:260: DownloadWarning: Downloading: http://naciscdn.org/naturalearth/110m/physical/ne_110m_coastline.zip warnings.warn('Downloading: {}'.format(url), DownloadWarning)

Exception occurred: File "/usr/share/miniconda/envs/xarray-tests/lib/python3.7/urllib/request.py", line 1319, in do_open raise URLError(err) urllib.error.URLError: <urlopen error [Errno 110] Connection timed out> The full traceback has been saved in /tmp/sphinx-err-nq73diee.log, if you want to report the issue to the developers. Please also report this if it was a user error, so that a better error message can be provided next time. A bug report can be filed in the tracker at https://github.com/sphinx-doc/sphinx/issues. Thanks!

[error]Bash exited with code '2'.

[section]Finishing: Build HTML docs

```

Dataset.copy(deep=True) and DataArray.copy (deep=True/False) accidentally cast IndexVariable's with dtype='<U*' to object. Same applies to copy.copy() and copy.deepcopy().

This is a regression in xarray >= 0.12.2. xarray 0.12.1 and earlier are unaffected.

```

In [1]: ds = xarray.Dataset( ...: coords={'x': ['foo'], 'y': ('x', ['bar'])}, ...: data_vars={'z': ('x', ['baz'])})

In [2]: ds
Out[2]: <xarray.Dataset> Dimensions: (x: 1) Coordinates: * x (x) <U3 'foo' y (x) <U3 'bar' Data variables: z (x) <U3 'baz'

In [3]: ds.copy()
Out[3]: <xarray.Dataset> Dimensions: (x: 1) Coordinates: * x (x) <U3 'foo' y (x) <U3 'bar' Data variables: z (x) <U3 'baz'

In [4]: ds.copy(deep=True)
Out[4]: <xarray.Dataset> Dimensions: (x: 1) Coordinates: * x (x) object 'foo' y (x) <U3 'bar' Data variables: z (x) <U3 'baz'

In [5]: ds.z
Out[5]: <xarray.DataArray 'z' (x: 1)> array(['baz'], dtype='<U3') Coordinates: * x (x) <U3 'foo' y (x) <U3 'bar'

In [6]: ds.z.copy()
Out[6]: <xarray.DataArray 'z' (x: 1)> array(['baz'], dtype='<U3') Coordinates: * x (x) object 'foo' y (x) <U3 'bar'

In [7]: ds.z.copy(deep=True)
Out[7]: <xarray.DataArray 'z' (x: 1)> array(['baz'], dtype='<U3') Coordinates: * x (x) object 'foo' y (x) <U3 'bar' ```

As of - dask 2.1.0 - distributed 2.1.0 - xarray 0.12.1 or git head (didn't try older versions):

```python

import xarray import distributed client = distributed.Client(set_as_default=False) ds = xarray.Dataset({'d': ('x', [1, 2])}).chunk(1) client.compute(ds).result() <xarray.Dataset> Dimensions: (x: 2) Dimensions without coordinates: x Data variables: d (x) int64 1 2

client.compute(ds.d).result() distributed.worker - WARNING - Compute Failed Function: _dask_finalize args: ([[array([1]), array([2])]], <function Dataset._dask_postcompute at 0x316a1db70>, ([(True, <this-array>, (<function Variable._dask_finalize at 0x3168f7f28>, (<function finalize at 0x1166bb8c8>, (), ('x',), OrderedDict(), None)))], set(), {'x': 2}, None, None, None, None), 'd') kwargs: {} Exception: KeyError(<this-array>)

KeyError Traceback (most recent call last) <ipython-input-8-2dbfe1b2ff17> in <module> ----> 1 client.compute(ds.d).result()

/anaconda3/lib/python3.7/site-packages/distributed/client.py in result(self, timeout) 226 result = self.client.sync(self._result, callback_timeout=timeout, raiseit=False) 227 if self.status == "error": --> 228 six.reraise(*result) 229 elif self.status == "cancelled": 230 raise result

/anaconda3/lib/python3.7/site-packages/six.py in reraise(tp, value, tb) 690 value = tp() 691 if value.traceback is not tb: --> 692 raise value.with_traceback(tb) 693 raise value 694 finally:

~/PycharmProjects/xarray/xarray/core/dataarray.py in _dask_finalize() 706 def _dask_finalize(results, func, args, name): 707 ds = func(results, *args) --> 708 variable = ds._variables.pop(_THIS_ARRAY) 709 coords = ds._variables 710 return DataArray(variable, coords, name=name, fastpath=True) ```

In variable.py, line 1116, you're missing a raise statement:

if isinstance(self.data, dask_array_type): TypeError("quantile does not work for arrays stored as dask " "arrays. Load the data via .compute() or .load() prior " "to calling this method.")

Currently looking into extending dask.percentile() to support more than 1D arrays, and then use it in xarray too.

Python 3.5.0 has gone out of the conda-forge repository. 3.5.1 is still there... for now. The anaconda repository starts directly from 3.5.4. 3.5.0 and 3.5.1 are a colossal pain in the back for typing support. Is this a good time to increase the requirement to >= 3.5.2? I honestly can't think how anybody could be unable to upgrade to the latest available 3.5 with minimal effort...

Can anybody tell me if there is anybody who benefits from this behaviour? I can't think of any good use cases.

``` wallet = xarray.DataArray([50, 70], dims=['currency'], coords={'currency': ['EUR', 'USD']}) restaurant_bill = xarray.DataArray([30], dims=['currency'], coords={'currency': ['USD']}) with xarray.set_options(arithmetic_join="outer"): print(wallet - restaurant_bill)

<xarray.DataArray (currency: 2)> array([ nan, 40.]) Coordinates: * currency (currency) object 'EUR' 'USD' ```

While it is fairly clear why it can be desirable to have nan + not nan = nan as a default in arithmetic when the nan is already present in one of the input arrays, when the nan is introduced as part of an automatic align things become much less intuitive.

Proposal: - add a parameter to xarray.align, fillvalue=numpy.nan, which determines what will appear in the newly created array elements - change __add__, __sub__ etc. to invoke xarray.align(fillvalue=0) - change __mul__, __truediv__ etc. to invoke xarray.align(fillvalue=1)

In theory the setting could be left as an opt-in as set_options(arithmetic_align_fillvalue='neutral'), yet I wonder who would actually want the current behaviour?

I'm using xarray to aggregate 38 GB worth of NetCDF data into a bunch of CSV reports. I have two problems:

1. The reports are 500,000 rows by 2,000 columns. Before somebody says "if you're using CSV for this size of data you're doing it wrong" - yes, I know, but it was the only way to make the data accessible to a bunch of people that only know how to use Excel and VBA. :tired_face: The sheer size of the reports means that (1) it's unsavory to keep the whole thing in RAM (2) pandas to_csv will take ages to complete (as it's single-threaded). The slowness is compounded by the fact that I have to compress everything with gzip.
2. I have to produce up to 40 reports from the exact same NetCDF files. I use dask to perform the computation, and different reports share a large amount of intermediate graph nodes. So I need to do everything in a single invocation to dask.compute() to allow the dask scheduler to de-duplicate the nodes.

To solve both problems, I wrote a new function: http://xarray-extras.readthedocs.io/en/latest/api/csv.html

And now my high level wrapper code looks like this: ```

DataSet from 200 .nc files, with a total of 500000 points on the 'row' dimension

nc = xarray.open_mfdataset('inputs..nc') reports = [ # DataArrays with shape (500000, 2000), with the rows split in 200 chunks gen_report0(nc), gen_report1(nc), .... gen_report39(nc), ] futures = [ # dask.delayed objects to_csv(reports[0], 'report0.csv.gz', compression='gzip'), to_csv(reports[1], 'report1.csv.gz', compression='gzip'), .... to_csv(reports[39], 'report39.csv.gz', compression='gzip'), ] dask.compute(futures) ``` The function is currently production quality in xarray-extras, but it would be very easy to refactor it as a method of xarray.DataArray in the main library.

Opinions?

DataArray.unstack() sorts the data alphabetically by label. Besides being poor for performance, this is very problematic whenever the order matters, and the labels are not in alphabetical order to begin with.

``` python

import xarray import pandas

index = [ ['x1', 'first' ], ['x1', 'second'], ['x1', 'third' ], ['x1', 'fourth'], ['x0', 'first' ], ['x0', 'second'], ['x0', 'third' ], ['x0', 'fourth'], ] index = pandas.MultiIndex.from_tuples(index, names=['x', 'count']) s = pandas.Series(list(range(8)), index) a = xarray.DataArray(s) a ```

<xarray.DataArray (dim_0: 8)> array([0, 1, 2, 3, 4, 5, 6, 7], dtype=int64) Coordinates: * dim_0 (dim_0) object ('x1', 'first') ('x1', 'second') ('x1', 'third') ...

python a.unstack('dim_0')

<xarray.DataArray (x: 2, count: 4)> array([[4, 7, 5, 6], [0, 3, 1, 2]], dtype=int64) Coordinates: * x (x) object 'x0' 'x1' * count (count) object 'first' 'fourth' 'second' 'third'

Whe the stack() method is used on a xarray with dask backend, one would expect that every output chunk is produced by exactly 1 input chunk.

This is not the case, as stack() actually produces an extremely fragmented dask array: https://gist.github.com/crusaderky/07991681d49117bfbef7a8870e3cba67

Is this by design? I can't put any sense in it ```

a = xarray.DataArray([1, 2, 3], dims=['x'], coords={'x': [1, 2, 3], 'y': 10}) a.coords['x'] <xarray.DataArray 'x' (x: 3)> array([1, 2, 3]) Coordinates: * x (x) int64 1 2 3 y int64 10 ```

In the xarray docs I read

With the current version of dask, there is no automatic alignment of chunks when performing operations between dask arrays with different chunk sizes. If your computation involves multiple dask arrays with different chunks, you may need to explicitly rechunk each array to ensure compatibility.

While chunk auto-alignment could be done within the dask library, that would be limited to arrays with the same dimensionality and same dims order. For example it would not be possible to have a dask library call to align the chunks on xarrays with the following dims: - (time, latitude, longitude) - (time) - (longitude, latitude)

even if it makes perfect sense in xarray.

I think xarray.align() should take care of it automatically.

A safe algorithm would be to always scale down the chunksize when in conflict. This would prevent having chunks larger than expected, and should minimise (in a greedy way) the number of operations. It's also a good idea on dask.distributed, where merging two chunks could cause one of them to travel on the network - which is very expensive.

e.g. to reconcile chunksizes a: (5, 10, 6) b: (5, 7, 9) the algorithm would rechunk both arrays to (5, 7, 3, 6).

Finally, when served with a numpy-based array and a dask-based array, align() should convert the numpy array to dask. The critical use case that would benefit from this behaviour is when align() is invoked inside a broadcast() between a tiny constant you just loaded from csv/pandas/pure python list/whatever - e.g. dims=(time, ) shape=(100, ) - and a huge dask-backed array e.g. dims=(time, scenario) shape=(100, 2**30) chunks=(25, 2**20).

Two unit tests are failing in the latest git master: - GenericNetCDFDataTest.test_cross_engine_read_write_netcdf3 - GenericNetCDFDataTestAutocloseTrue.test_cross_engine_read_write_netcdf3

Both with the message:

``` xarray/tests/test_backends.py:1558:

xarray/backends/api.py:286: in open_dataset autoclose=autoclose) xarray/backends/netCDF4_.py:275: in open ds = opener() xarray/backends/netCDF4_.py:199: in _open_netcdf4_group ds = nc4.Dataset(filename, mode=mode, **kwargs) netCDF4/_netCDF4.pyx:2015: in netCDF4._netCDF4.Dataset.init ???

??? E OSError: [Errno -36] NetCDF: Invalid argument: b'/tmp/tmpwp675lnc/temp-1069.nc'

netCDF4/_netCDF4.pyx:1636: OSError ```

Attaching conda list: conda.txt

```

ds = xarray.Dataset({'x': 1}) ds.to_netcdf('foo.nc') dask/utils.py:1010: UserWarning: Deprecated, see dask.base.get_scheduler instead ```

Stack trace: ```

xarray/backends/common.py(44)get_scheduler() 43 from dask.utils import effective_get ---> 44 actual_get = effective_get(get, collection) ``` There are two separate problems here:

• dask recently changed API from get(get=callable) to get(scheduler=str). Should we
• just increase the minimum version of dask (I doubt anybody will complain)
• go through the hoops of dynamically invoking a different API depending on the dask version :sweat:
• silence the warning now, and then increase the minimum version of dask the day that dask removes the old API entirely (risky)?
• xarray is calling dask even when it's unnecessary, as none of the variables in the example Dataset had a dask backend. I don't think there are any CI suites for NetCDF without dask. I'm also wondering if they would bring any actual added value, as dask is small, has no exotic dependencies, and is pure Python; so I doubt anybody will have problems installing it whatever his setup is.

@shoyer opinion?

In xarray 0.10.4, the dtype encoding in to_netcdf has stopped working, for all engines: ```

import xarray ds = xarray.Dataset({'x': ['foo', 'bar', 'baz']}) ds.to_netcdf('test.nc', encoding={'x': {'dtype': 'S1'}}) [...]

xarray/backends/netCDF4_.py in _extract_nc4_variable_encoding(variable, raise_on_invalid, lsd_okay, h5py_okay, backend, unlimited_dims) 196 if invalid: 197 raise ValueError('unexpected encoding parameters for %r backend: ' --> 198 ' %r' % (backend, invalid)) 199 else: 200 for k in list(encoding):

ValueError: unexpected encoding parameters for 'netCDF4' backend: ['dtype'] ``` I'm still trying to figure out how the regression tests didn't pick it up and what change introduced it.

@shoyer I'm working on this as my top priority. Do you agree this is serious enough for an emergency re-release? (0.10.4.1 or 0.10.5, your choice)

Dataset().update(Dataset())

FutureWarning: iteration over an xarray.Dataset will change in xarray v0.11 to only include data variables, not coordinates. Iterate over the Dataset.variables property instead to preserve existing behavior in a forwards compatible manner.

This is a regression in xarray 0.10.4. @shoyer this isn't serious enough to warrant an immediate release on its own, but we're already doing one so we might as well include it.

In this script: ``` import pandas import xarray

df = pandas.DataFrame( [[1, 2], [3, 4]], index=['foo', 'foo'], columns=['bar', 'baz']) print(df.stack())

a = xarray.DataArray(df) print(a.stack(s=a.dims)) ```

The first part works both with pandas 0.22 and 0.23. The second part works in xarray 0.10. 4 + pandas 0.22, and crashes with pandas 0.23:

File "/mnt/resource/tmp/anaconda_guido/lib/python3.6/site-packages/xarray/core/dataarray.py", line 1115, in stack ds = self._to_temp_dataset().stack(**dimensions) File "/mnt/resource/tmp/anaconda_guido/lib/python3.6/site-packages/xarray/core/dataset.py", line 2123, in stack result = result._stack_once(dims, new_dim) File "/mnt/resource/tmp/anaconda_guido/lib/python3.6/site-packages/xarray/core/dataset.py", line 2092, in _stack_once idx = utils.multiindex_from_product_levels(levels, names=dims) File "/mnt/resource/tmp/anaconda_guido/lib/python3.6/site-packages/xarray/core/utils.py", line 96, in multiindex_from_product_levels return pd.MultiIndex(levels, labels, sortorder=0, names=names) File "/mnt/resource/tmp/anaconda_guido/lib/python3.6/site-packages/pandas/core/indexes/multi.py", line 240, in __new__ result._verify_integrity() File "/mnt/resource/tmp/anaconda_guido/lib/python3.6/site-packages/pandas/core/indexes/multi.py", line 283, in _verify_integrity level=i)) ValueError: Level values must be unique: ['foo', 'foo'] on level 0

As of today, Dataset.to_netcdf() exclusively allows writing uncompressed or compressed with zlib. zlib was absolutely revolutionary when it was released... in 1995. Time has passed, and much better compression algorithms have appeared over time. Good news is, h5py supports LZF out of the box, and is extensible with plugins to support theoretically any other algorithm. h5netcdf exposes such interface through its new (non-legacy) API; however Dataset.to_netcdf(engine='h5netcdf') supports the legacy API exclusively.

I already tested that, once you manage to write to disk with LZF (using h5netcdf directly), open_dataset(engine='h5netcdf') transparently opens the compressed store.

Options: - write a new engine for Dataset.to_netcdf() to support the new h5netcdf API. - switch the whole engine='h5netcdf' to the new API. Drop support for the old parameters in to_netcdf(). This is less bad than it sounds, as people can switch to another engine in case of trouble. This is the cleanest solution, but also the most disruptive one. - switch the whole engine='h5netcdf' to the new API; have to_netcdf() accept both new and legacy parameters, and implement a translation layer of parameters from the legacy API to the new API. The benefit here is that, as long as the user sticks to the legacy API, he can hop between engines transparently. On the other hand I have a hard time believing anybody would care. - ?

I've written a series of wrappers for the 1-dimensional scipy interpolators.

Prototype code and colourful demo plots: https://gist.github.com/crusaderky/b0aa6b8fdf6e036cb364f6f40476cc67

Features

• Interpolate a ND array on any arbitrary dimension
• Nearest-neighbour, linear, quadratic, cubic, Akima, PCHIP, and custom interpolators are supported
• dask supported on both on the interpolated array and x_new
• Supports ND x_new arrays
• The CPU-heavy interpolator generation (splrep) is executed only once and then can be applied to multiple x_new (splev)
• Pickleable and distributed-friendly

Design hacks

• Depends on dask module, even when all inputs are based on plain numpy.
• Abuses attrs and the ability to invoke a.attrname to get the user experience of a new DataArray method.
• Abuses the fact that the chunks of a dask.array.Array can contain anything and you won't notice until you compute them.

Limitations

• Can't dump to netcdf. Not solvable without hacking into the implementation details of scipy.
• Datasets are not supported. Trivial to fix after solving #1699.
• Chunks are not supported on x_new. Trivial to fix after solving #1995.
• Chunks are not supported along the interpolated dimension. This is very complicated to solve. If x and x_new were always monotonic ascending,it would be (not trivially) solvable with dask.array.ghost.ghost. If you make no assumptions about monotonicity, things become way more complicated. A solution would need to go in the dask module, and then be invoked trivially from here with dask='allowed'.

xarray.dot() has comparable performance with numpy.einsum. However, when it uses a dask backend, it's much slower than the new dask.array.einsum function (https://github.com/dask/dask/pull/3412). The performance gap widens when the dimension upon which you are reducing is chunked.

Also, for some reason dot(a<s, t>, b<t>, dims=[t]) and dot(a<s,t>, a<s,t>, dims=[s,t]) do work (very slowly) when s and t are chunked, while dot(a<s, t>, a<s, t>, dims=[t]) crashes complaining it can't operate on a chunked core dim (related discussion: https://github.com/pydata/xarray/issues/1995).

The proposed solution is to simply wait for https://github.com/dask/dask/pull/3412 to reach the next release and then reimplement xarray.dot to use dask.array.einsum. This means that dask users will lose the ability to use xarray.dot if they upgrade xarray version but not dask version, but I believe it shouldn't be a big problem for most?

``` import numpy import dask.array import xarray

def bench(tchunk, a_by_a, dims, iis): print(f"\nbench({tchunk}, {a_by_a}, {dims}, {iis})")

a = xarray.DataArray(
    dask.array.random.random((500000, 100), chunks=(50000, tchunk)),
    dims=['s', 't'])
if a_by_a:
    b = a
else:
    b = xarray.DataArray(
        dask.array.random.random((100, ), chunks=tchunk),
        dims=['t'])

print("xarray.dot(numpy backend):")
%timeit xarray.dot(a.compute(), b.compute(), dims=dims)
print("numpy.einsum:")
%timeit numpy.einsum(iis, a, b)
print("xarray.dot(dask backend):")
try:
    %timeit xarray.dot(a, b, dims=dims).compute()
except ValueError as e:
    print(e)
print("dask.array.einsum:")
%timeit dask.array.einsum(iis, a, b).compute()

bench(100, False, ['t'], '...i,...i') bench( 20, False, ['t'], '...i,...i') bench(100, True, ['t'], '...i,...i') bench( 20, True, ['t'], '...i,...i') bench(100, True, ['s', 't'], '...ij,...ij') bench( 20, True, ['s', 't'], '...ij,...ij') Output: bench(100, False, ['t'], ...i,...i) xarray.dot(numpy backend): 195 ms ± 3.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) numpy.einsum: 205 ms ± 2.47 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) xarray.dot(dask backend): 356 ms ± 44.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) dask.array.einsum: 244 ms ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

bench(20, False, ['t'], ...i,...i) xarray.dot(numpy backend): 297 ms ± 16.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) numpy.einsum: 254 ms ± 15.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) xarray.dot(dask backend): 732 ms ± 74.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) dask.array.einsum: 274 ms ± 12.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

bench(100, True, ['t'], ...i,...i) xarray.dot(numpy backend): 438 ms ± 43.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) numpy.einsum: 415 ms ± 17.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) xarray.dot(dask backend): 633 ms ± 31.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) dask.array.einsum: 431 ms ± 17 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

bench(20, True, ['t'], ...i,...i) xarray.dot(numpy backend): 457 ms ± 17.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) numpy.einsum: 463 ms ± 24.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) xarray.dot(dask backend): dimension 't' on 0th function argument to apply_ufunc with dask='parallelized' consists of multiple chunks, but is also a core dimension. To fix, rechunk into a single dask array chunk along this dimension, i.e., .rechunk({'t': -1}), but beware that this may significantly increase memory usage. dask.array.einsum: 485 ms ± 15.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

bench(100, True, ['s', 't'], ...ij,...ij) xarray.dot(numpy backend): 418 ms ± 14.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) numpy.einsum: 444 ms ± 43.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) xarray.dot(dask backend): 384 ms ± 57.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) dask.array.einsum: 415 ms ± 19.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

bench(20, True, ['s', 't'], ...ij,...ij) xarray.dot(numpy backend): 489 ms ± 2.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) numpy.einsum: 443 ms ± 3.35 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) xarray.dot(dask backend): 585 ms ± 64.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) dask.array.einsum: 455 ms ± 13.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) ```

Algorithms that are wrapped by xarray.apply_ufunc(dask='parallelized'), and in general most algorithms for which aren't embarassingly parallel and for which there isn't a sophisticated dask function that allows for multiple chunks, cannot have multiple chunks on their core dimensions.

I have lost count of how many times I prefixed my invocations of apply_ufunc on a DataArray with the same blurb, over and over again: if x.chunks: x = x.chunk({dim: x.shape[x.dims.index(dim)]}) The reason why it looks so awful is that DataArray.shape, DataArray.dims, Variable.shape and Variable.dims are positional.

I can see a few possible solutions to the problem:

Design 1

Change DataArray.chunk etc. to accept a special chunk size, e.g. -1, which means "whatever the size of that dim is". The above would become: if x.chunks: x = x.chunk({dim: -1}) which is much more bearable. One could argue that the implementation would need to happen in dask.array.rechunk; on the other hand in dask it woulf feel silly, because already today you can do it in a very synthetic way: x = x.rechunk({axis: x.shape[axis]}) I'm not overly fond of this solution as it would be rather obscure for anybody who isn't super familiar with the API documentation.

Design 2

Add properties to DataArray and Variable, ddims and dshape (happy to hear suggestions about better names), which would return dims and shape as a OrderedDict, just like Dataset.dims and Dataset.shape.

The above would become: if x.chunks: x = x.chunk({dim: x.dshape[dim]})

Design 3

Change dask.array.rechunk to accept numpy.inf / math.inf as the chunk size. This makes sense, as the function already accepts chunk sizes that are larger than the shape - however, it's currently limited to int. This is probably my personal favourite, and trivial to implement too.

The above would become: if x.chunks: x = x.chunk({dim: np.inf})

Design 4

Introduce a convenience method for DataArray, Dataset, and Variable, ensure_single_chunk(*dims). Below a prototype: ``` def ensure_single_chunk(a, *dims): """If a has dask backend and two or more chunks on dims, rechunk it so that they become single-chunked. This is typically a prerequisite for computing any algorithm along dim that is not embarassingly parallel (short of sophisticated implementations such as those found in the dask module).

:param a:
    any xarray object
:param str dims:
    one or more dims of a to rechunk
:returns:
    copy of a, where all listed dims are guaranteed to be on a single dask chunk.
    if a has numpy backend, return a shallow copy of it.
"""
if isinstance(a, xarray.Dataset):
    dims = set(dims)
    unknown_dims = dims - a.dims.keys()
    if unknown_dims:
        raise ValueError("dim(s) %s not found" % ",".join(unknown_dims))
    a = a.copy(deep=False)
    for k, v in a.variables.items():
        if v.chunks:
            a[k] = ensure_single_chunk(v, *(set(v.dims) & dims))
    return a

if not isinstance(a, (xarray.DataArray, xarray.Variable)):
    raise TypeError('a must be a DataArray, Dataset, or Variable')

if not a.chunks:
    # numpy backend
    return a.copy(deep=False)

return a.chunk({
    dim: a.shape[a.dims.index(dim)]
    for dim in dims
})

```

As of xarray-0.10-rc1:

Works: ``` import xarray import scipy.stats a = xarray.DataArray([1,2], dims=['x'])

xarray.apply_ufunc(scipy.stats.norm.cdf, a, 0, 1)

<xarray.DataArray (x: 2)> array([ 0.841345, 0.97725 ]) Dimensions without coordinates: x ```

Broken: ``` xarray.apply_ufunc( scipy.stats.norm.cdf, a.chunk(), 0, 1, dask='parallelized', output_dtypes=[a.dtype] ).compute()

IndexError Traceback (most recent call last) <ipython-input-35-1d4025e1ebdb> in <module>() ----> 1 xarray.apply_ufunc(scipy.stats.norm.cdf, a.chunk(), 0, 1, dask='parallelized', output_dtypes=[a.dtype]).compute()

~/anaconda3/lib/python3.6/site-packages/xarray/core/computation.py in apply_ufunc(func, args, kwargs) 913 join=join, 914 exclude_dims=exclude_dims, --> 915 keep_attrs=keep_attrs) 916 elif any(isinstance(a, Variable) for a in args): 917 return variables_ufunc(args)

~/anaconda3/lib/python3.6/site-packages/xarray/core/computation.py in apply_dataarray_ufunc(func, args, kwargs) 210 211 data_vars = [getattr(a, 'variable', a) for a in args] --> 212 result_var = func(data_vars) 213 214 if signature.num_outputs > 1:

~/anaconda3/lib/python3.6/site-packages/xarray/core/computation.py in apply_variable_ufunc(func, args, kwargs) 561 raise ValueError('unknown setting for dask array handling in ' 562 'apply_ufunc: {}'.format(dask)) --> 563 result_data = func(input_data) 564 565 if signature.num_outputs > 1:

~/anaconda3/lib/python3.6/site-packages/xarray/core/computation.py in <lambda>(arrays) 555 func = lambda arrays: _apply_with_dask_atop( 556 numpy_func, arrays, input_dims, output_dims, signature, --> 557 output_dtypes, output_sizes) 558 elif dask == 'allowed': 559 pass

~/anaconda3/lib/python3.6/site-packages/xarray/core/computation.py in _apply_with_dask_atop(func, args, input_dims, output_dims, signature, output_dtypes, output_sizes) 624 for element in (arg, dims[-getattr(arg, 'ndim', 0):])] 625 return da.atop(func, out_ind, *atop_args, dtype=dtype, concatenate=True, --> 626 new_axes=output_sizes) 627 628

~/anaconda3/lib/python3.6/site-packages/dask/array/core.py in atop(func, out_ind, args, kwargs) 2231 raise ValueError("Must specify dtype of output array") 2232 -> 2233 chunkss, arrays = unify_chunks(args) 2234 for k, v in new_axes.items(): 2235 chunkss[k] = (v,)

~/anaconda3/lib/python3.6/site-packages/dask/array/core.py in unify_chunks(args, *kwargs) 2117 chunks = tuple(chunkss[j] if a.shape[n] > 1 else a.shape[n] 2118 if not np.isnan(sum(chunkss[j])) else None -> 2119 for n, j in enumerate(i)) 2120 if chunks != a.chunks and all(a.chunks): 2121 arrays.append(a.rechunk(chunks))

~/anaconda3/lib/python3.6/site-packages/dask/array/core.py in <genexpr>(.0) 2117 chunks = tuple(chunkss[j] if a.shape[n] > 1 else a.shape[n] 2118 if not np.isnan(sum(chunkss[j])) else None -> 2119 for n, j in enumerate(i)) 2120 if chunks != a.chunks and all(a.chunks): 2121 arrays.append(a.rechunk(chunks))

IndexError: tuple index out of range ```

Workaround: ``` xarray.apply_ufunc( scipy.stats.norm.cdf, a, kwargs={'loc': 0, 'scale': 1}, dask='parallelized', output_dtypes=[a.dtype]).compute()

<xarray.DataArray (x: 2)> array([ 0.841345, 0.97725 ]) Dimensions without coordinates: x ```

I have 200x of the below dataset, split on the 'scenario' axis:

<xarray.Dataset> Dimensions: (fx_id: 39, instr_id: 16095, scenario: 2501) Coordinates: currency (instr_id) object 'GBP' 'USD' 'GBP' 'GBP' 'GBP' 'EUR' 'CHF' ... * fx_id (fx_id) object 'USD' 'EUR' 'JPY' 'ARS' 'AUD' 'BRL' 'CAD' ... * instr_id (instr_id) object 'property_standard_gbp' ... * scenario (scenario) object 'Base Scenario' 'SSMC_1' 'SSMC_2' ... type (instr_id) object 'Common Stock' 'Fixed Amortizing Bond' ... Data variables: fx_rates (fx_id, scenario) float64 1.236 1.191 1.481 1.12 1.264 ... instruments (instr_id, scenario) float64 1.0 1.143 0.9443 1.013 1.176 ... Attributes: base_currency: GBP

I individually dump them to disk with Dataset.to_netcdf(fname, engine='h5netcdf'). Then I try loading them back up with open_mfdataset, but it's mortally slow:

``` %%time xarray.open_mfdataset('*.nc', engine='h5netcdf')

Wall time: 30.3 s ```

The problem is caused by the coords being read from disk multiple times. Workaround:

%%time def load_coords(ds): for coord in ds.coords.values(): coord.load() return ds xarray.open_mfdataset('*.nc', engine='h5netcdf', preprocess=load_coords) Wall time: 12.3 s

Proposed solutions: 1. Implement the above workaround directly inside open_mfdataset() 2. change open_dataset() to always eagerly load the coords to memory, regardless of the chunks parameter. Is there any valid use case where lazy coords are actually desirable?

An additional, more radical observation is that, very frequently, a user knows in advance that all coords are aligned. In this use case, the user could explicitly request xarray to blindly trust this assumption, and thus skip loading the coords not based on concat_dim in all datasets beyond the first.

ds = Dataset() ds.encoding = {"unlimited_dims": 'x'} ds.copy().encoding {} By looking at dataset.py, there's a lot of calls to Dataset._construct_direct that omit the encoding. Is it correct to add it in all cases?

I think I'm losing my sanity here. I have a anaconda3 Python 3.6 environment with all required and optional dependencies of xarray installed and updated to the latest available version, except pyNio. If I run test.py on the latest xarray package from the git tip, the vast majority of the tests in test_backends.py are skipped - including those that have nothing to do with pyNio! e.g.

tests/test_backends.py::ScipyInMemoryDataTest::test_bytesio_pickle PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_coordinates_encoding SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_dataset_caching SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_dataset_compute SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_default_fill_value SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_encoding_kwarg SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_encoding_same_dtype SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_invalid_dataarray_names_raise SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_load SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_orthogonal_indexing PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_pickle SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_pickle_dataarray SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_None_variable SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_boolean_dtype SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_coordinates SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_datetime_data SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_endian SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_example_1_netcdf SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_float64_data SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_mask_and_scale SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_object_dtype SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_string_data SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_strings_with_fill_value SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_test_data SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_timedelta_data SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_to_netcdf_explicit_engine PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_unsigned_roundtrip_mask_and_scale SKIPPED tests/test_backends.py::ScipyInMemoryDataTest::test_write_store PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_zero_dimensional_variable SKIPPED

If I comment out line 1462:

``` @requires_scipy

@requires_pynio

class TestPyNio(CFEncodedDataTest, Only32BitTypes, TestCase): ```

Then magically everything starts working again!

tests/test_backends.py::ScipyInMemoryDataTest::test_bytesio_pickle PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_coordinates_encoding PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_dataset_caching PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_dataset_compute PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_default_fill_value PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_encoding_kwarg PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_encoding_same_dtype PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_invalid_dataarray_names_raise PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_load PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_orthogonal_indexing PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_pickle PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_pickle_dataarray PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_None_variable PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_boolean_dtype PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_coordinates PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_datetime_data PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_endian PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_example_1_netcdf PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_float64_data PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_mask_and_scale PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_object_dtype PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_string_data PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_strings_with_fill_value PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_test_data PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_roundtrip_timedelta_data PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_to_netcdf_explicit_engine PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_unsigned_roundtrip_mask_and_scale PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_write_store PASSED tests/test_backends.py::ScipyInMemoryDataTest::test_zero_dimensional_variable PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_bytesio_pickle PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_coordinates_encoding PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_dataset_caching PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_dataset_compute PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_default_fill_value PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_encoding_kwarg PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_encoding_same_dtype PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_invalid_dataarray_names_raise PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_load PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_orthogonal_indexing PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_pickle PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_pickle_dataarray PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_None_variable PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_boolean_dtype PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_coordinates PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_datetime_data PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_endian PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_example_1_netcdf PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_float64_data PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_mask_and_scale PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_object_dtype PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_string_data PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_strings_with_fill_value PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_test_data PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_roundtrip_timedelta_data PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_to_netcdf_explicit_engine PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_unsigned_roundtrip_mask_and_scale PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_write_store PASSED tests/test_backends.py::ScipyInMemoryDataTestAutocloseTrue::test_zero_dimensional_variable PASSED

duck_array_ops.concatenate and duck_array_ops.stack load dask variables if the first one is numpy-based:

``` xarray.concat([ xarray.DataArray([1]).chunk(), xarray.DataArray([1]), ], dim='dim_0')

Out[1]: <xarray.DataArray (dim_0: 2)> dask.array<shape=(2,), dtype=int64, chunksize=(1,)> Dimensions without coordinates: dim_0

xarray.concat([ xarray.DataArray([1]), xarray.DataArray([1]).chunk(), ], dim='dim_0')

Out[2]: <xarray.DataArray (dim_0: 2)> array([1, 1]) Dimensions without coordinates: dim_0 ```

DataArray.__repr__ and Variable.__repr__ print a placeholder if the data uses the dask backend. Not so Dataset.__repr__, which tries computing the data before printing a tiny preview of it. This issue is extremely annoying when working in Jupyter, and particularly acute if the chunks are very big or are at the end of a very long chain of computation.

For data variables, the expected behaviour is to print a placeholder just like DataArray does. For coords, we could either - print a placeholders (same treatment as data variables) - automatically invoke load() when the coord is added to the dataset - see #1521 for discussion.

I work with a very large dask-based algorithm in xarray, and I do my optimization by hand before hitting compute(). In other cases, I need using multiple dask schedulers at once (e.g. a multithreaded one for numpy-based work and a multiprocessing one for pure python work).

This change proposal (which I'm happy to do) is about accepting *args, **kwds parameters in all .compute(), .load(), and .persist() xarray methods and pass them verbatim to the underlying dask compute() and persist() functions.

Whenever xarray performs a shallow copy of any object (DataArray, Dataset, Variable), it creates a view of the underlying numpy arrays. This design fails when the object is pickled.

Whenever a numpy view is pickled, it becomes a regular array:

```

a = numpy.arange(226) print(len(pickle.dumps(a)) / 220) 256.00015354156494 b = a.view() print(len(pickle.dumps((a, b))) / 2**20) 512.0001964569092 b.base is a True a2, b2 = pickle.loads(pickle.dumps((a, b))) b2.base is a2 False ```

This has devastating effects in my use case. I start from a dask-backed DataArray with a dimension of 500,000 elements and no coord, so the coord is auto-assigned by xarray as an incremental integer. Then, I perform ~3000 transformations and dump the resulting dask-backed array with pickle. However, I have to dump all intermediate steps for audit purposes as well. This means that xarray invokes numpy.arange to create (500k * 4 bytes) ~ 2MB worth of coord, then creates 3000 views of it, which the moment they're pickled expand to several GBs as they become 3000 independent copies.

I see a few possible solutions to this: 1. Implement pandas range indexes in xarray. This would be nice as a general thing and would solve my specific problem, but anybody who does not fall in my very specific use case won't benefit from it. 2. Do not auto-generate a coord with numpy.arange() if the user doesn't explicitly ask for it; just leave a None and maybe generate it on the fly when requested. Again, this would solve my specific problem but not other people's. 3. Force the coord to be a dask.array.arange. Actually supporting unconverted dask arrays as coordinates would take a considerable amount of work; they would get converted to numpy several times, and other issues. Again it wouldn't solve the general problem. 4. Fix the issue upstream in numpy. I didn't look into it yet and it's definitely worth investigating, but I found about it as early as 2012, so I suspect there might be some pretty good reason why it works like that... 5. Whenever xarray performs a shallow copy, take the numpy array instead of creating a view.

I implemented (5) as a workaround in my getstate method. Before:

%%time print(len(pickle.dumps(cache, pickle.HIGHEST_PROTOCOL)) / 2**30) 2.535497265867889 Wall time: 33.3 s

Workaround:

``` def get_base(array): if not isinstance(array, numpy.ndarray): return array
elif array.base is None: return array elif array.base.dtype != array.dtype: return array elif array.base.shape != array.shape: return array else: return array.base

for v in cache.values(): if isinstance(v, xarray.DataArray): v.data = get_base(v.data) for coord in v.coords.values(): coord.data = get_base(coord.data) elif isinstance(v, xarray.Dataset): for var in v.variables(): var.data = get_base(var.data) ```

After:

%%time print(len(pickle.dumps(cache, pickle.HIGHEST_PROTOCOL)) / 2**30) 0.9733252348378301 Wall time: 21.1 s

``` python

a = xarray.DataArray([1,2]).chunk(1) a <xarray.DataArray (dim_0: 2)> dask.array<xarray-..., shape=(2,), dtype=int64, chunksize=(1,)> Coordinates: * dim_0 (dim_0) int64 0 1 xarray.broadcast(a) (<xarray.DataArray (dim_0: 2)> array([1, 2]) Coordinates: * dim_0 (dim_0) int64 0 1,) ```

The problem is actually somewhere in the constructor of DataArray. In alignment.py:362, we have return DataArray(data, ...) where data is a Variable with dask backend. The returned DataArray object has a numpy backend. As a workaround, changing that line to return DataArray(data.data, ...) (thus passing a dask array) fixes the problem.

After that however there's a new issue: whenever broadcast adds a dimension to an array, it creates it in a single chunk, as opposed to copying the chunking of the other arrays. This can easily call a host to go out of memory, and makes it harder to work with the arrays afterwards because chunks won't match.

I'd like to add the following top-level functions to xarray:

``` def const_like(array, value=0): """Return a new array with the same shape of array and the given constant value. If array is dask-backed, return a new dask-backed array with the same chunks.

:param array:
    a numpy or dask-backed xarray.DataArray
:param value:
    any scalar number
"""
if isinstance(array.data, dask.array.Array):
    if value == 0:
        data = dask.array.zeros(
            array.data.shape,
            chunks=array.data.chunks,
            dtype=array.data.dtype)
    else:
        data = dask.array.ones(
            array.data.shape,
            chunks=array.data.chunks,
            dtype=array.data.dtype)
else:
    if value == 0:
        data = numpy.zeros_like(array.data)
    else:
        data = numpy.ones_like(array.data)
if value not in (0, 1):
    data = data * value

return xarray.DataArray(data, dims=array.dims, coords=array.coords, attrs=array.attrs)

def zeros_like(array): return const_like(array, 0)

def ones_like(array): return const_like(array, 1) ```

The above would need to be expanded to support Dataset and Variable objects. In Datasets, the data_vars would be constants whereas all other variables would be copied verbatim. Thoughts?

Pickling a xarray.DataArray with dask backend will cause it to resolve the .data to a numpy array. This is not desirable, as there's legitimate use cases where you may want to e.g. save a computation for later, or send it somewhere across the network.

Analogously, auto-converting a dask xarray to a numpy xarray as soon as you invoke the .value property is probably nice when you are working on a jupyter terminal, but not in a general purpose situation, particularly when xarray is used at the foundation of a very complex framework. Most of my headaches so far have been caused trying to figure out when, where and why the dask backend was replaced with numpy.

IMHO a module-wide switch to disable implicit dask->numpy conversion would be a nice solution. A new method, compute(), could explicitly convert in place from dask to numpy.

I have two arrays with misaligned dimensions x and y, and I want to concatenate them on dimension y. I can't seem to find any way to do it, because: 1. If I do not invoke align(), it will fail complaining that dimension x is not aligned 2. if I invoke align(), it will create unwanted elements on dimension y

See example: https://gist.github.com/crusaderky/a96db5b59396d94fe1e22694bc091d55

Am I missing something obvious? Possibly align() should accept an optional parameter e.g. exclude==['y']?

Thanks in advance

It's currently impossible to pickle xarray.ufuncs.

import xarray.ufuncs, pickle pickle.dumps(xarray.ufuncs.maximum)

AttributeError: Can't pickle local object '_create_op.<locals>.func'

Take a dask-backed array: a = xarray.DataArray(dask.array.random.random(100 * 2**30, chunks=2**20)) This works: b = xarray.ufuncs.maximum(a, 0) This will cripple your computer and force you to reboot: b = xarray.ufuncs.maximum(0, a)

In the second case, xarray.ufuncs.maximum is resolving the dask array - in other wods, it's doing numpy.maximum(0, a.values)