issue_comments

But for this application in particular presumably all of these already tasks fuse together into a single composite task, yes? Then we would need to commute certain operations in some order (preferring to push down slicing operations). Then we would need to merge all of the slices, yes?

Yes, that sounds right to me.

Ah, so here's the thing: decode_cf is actually lazy, but it uses some old xarray machinery for lazily decoding conventions instead of dask: In [21]: xr.decode_cf(chunked).foo.variable._data Out[21]: LazilyIndexedArray(array=dask.array<xarray-foo, shape=(100,), dtype=int64, chunksize=(100,)>, key=(slice(None, None, None),)) You could call .chunk() on it again to load it into dask, but even though the whole thing is lazy, the new dask arrays don't know they are holding dask arrays, which makes life difficult for dask.

You might ask why this separate lazy compute machinery exists. The answer is that dask fails to optimize element-wise operations like (scale * array)[subset] -> scale * array[subset], which is a critical optimization for lazy decoding of large datasets.

See https://github.com/dask/dask/issues/746 for discussion and links to PRs about this. @jcrist had a solution that worked, but it slowed down every dask array operations by 20%, which wasn't a great win.

I wonder if this is worth revisiting with a simpler, less general optimization pass that doesn't bother with broadcasting. See the subclasses of NDArrayMixin in xarray/conventions.py for examples of the sorts of functionality we need: - Casting (e.g., array.astype(bool)). - Chained arithmetic with scalars (e.g., 0.5 + 0.5 * array). - Custom element-wise operations (e.g., map_blocks(convert_to_datetime64, array, dtype=np.datetime64)) - Custom aggregations that drop a dimension (e.g., map_blocks(characters_to_string, array, drop_axis=-1))

If we could optimize all these operations (and ideally chain them), then we could drop all the lazy loading stuff from xarray in favor of dask, which would be a real win. @mrocklin any thoughts?

A simple test case: ``` In [16]: ds = xr.Dataset({'foo': ('x', np.arange(100))})

In [17]: chunked = ds.chunk()

In [18]: chunked.foo Out[18]: <xarray.DataArray 'foo' (x: 100)> dask.array<xarray-foo, shape=(100,), dtype=int64, chunksize=(100,)> Dimensions without coordinates: x

In [19]: xr.decode_cf(chunked).foo Out[19]: <xarray.DataArray 'foo' (x: 100)> array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99]) Dimensions without coordinates: x ```

decode_cf() should definitely work lazily on dask array. If not, I would consider that a bug.