Ah ha! Excellent. Thanks @shoyer . I'll give this a shot tomorrow (or perhaps ask @jcrist to look into it if he has time).

As expected, the following all dask.array solution triggers this:

``` python dates = pd.date_range('2001-01-01', freq='D', periods=1000) sizes = pd.Series(dates, dates).resample('1M', how='count').values chunks = (tuple(sizes), (100,)) x = da.ones((3630, 100), chunks=chunks) assert x[240:270].shape == x[240:270].compute().shape

AssertionError

```

Something like this might work to generate pathological chunks for dask.array:

dates = pandas.date_range('2000-01-01', freq='D', periods=1000) sizes = pandas.Series(dates, dates).resample('1M', how='count').values chunks = (tuple(sizes), (100,))

(I don't have xarray or dask installed on my work computer, but I could check this later)

@shoyer perhaps you can help to translate the code within @pwolfram 's script (in particular the lines that I've highlighted) and say how xarray would use dask.array to accomplish this. rnum = 7, Ntr = 30

I think this is a case where we each have some necessary expertise to resolve this issue. We probably need to work together to efficiently hunt down what's going on.

If you don't specify a chunksize, xarray should use each file as a full "chunk". So it would probably be useful to know the shapes of each array you are loading with open_mfdataset. My guess is that this issue only arises when indexing arrays consisting of differently sized chunks, which is exactly why using .chunk to set a fixed chunk size resolves this issue.

To be clear, all the logic implementing the chunking and indexing code for xarray objects containing dask arrays lives inside dask.array itself, not in our xarray wrapper (which is pretty thin). This doesn't make this any less of an issue for you, but I'm pretty sure (and I think @mrocklin agrees) that the bug here in probably in the dask.array layer.

It looks like the issue is in these lines:

(Pdb) pp rlzns.xParticle.data dask.array<getitem..., shape=(3630, 100), dtype=float64, chunksize=(21, 100)> (Pdb) pp rlzns.xParticle[rnum*Ntr:(rnum+1)*Ntr,:].data dask.array<getitem..., shape=(30, 100), dtype=float64, chunksize=(23, 100)> (Pdb) pp rlzns.xParticle[rnum*Ntr:(rnum+1)*Ntr,:].data.compute().shape (29, 100)

I'm confused by the chunksize change from 21 to 23.

In straight dask.array I'm unable to reproduce this problem, although obviously I'm doing something differently here than how xarray does things.

``` python In [1]: import dask.array as da x In [2]: x = da.ones((3630, 100), chunks=(21, 100))

In [3]: y = x[730:830, :]

In [4]: y.shape Out[4]: (30, 100)

In [5]: y.compute().shape Out[5]: (30, 100)

In [6]: y.chunks Out[6]: ((21, 9), (100,)) ```

It would be awesome if you all could produce a failing example with just dask.array.

This does look very strange. I'm guessing it's a dask.array bug (cc @mrocklin).

Can you make a reproducible example? If so, we'll probably be able to figure this out. How do you make this data?

Tracking this sort of thing down is a good motivation for an eager-evaluation mode in dask.array... (https://github.com/dask/dask/issues/292)