I have a large DataSet, including these DataArrays:

<xarray.DataArray 'temp' (ocean_time: 1325, s_rho: 30, eta_rho: 602, xi_rho: 677)>
dask.array<concatenate, shape=(1325, 30, 602, 677), dtype=float32, chunksize=(1, 1, 602, 677)>

and

<xarray.DataArray 'zeta' (ocean_time: 1325, eta_rho: 602, xi_rho: 677)>
dask.array<concatenate, shape=(1325, 602, 677), dtype=float32, chunksize=(1, 602, 677)>

(The coordinates and attributes excluded for brevity, but they match in the right ways.)

When I do math operations with the 4D DataArray (temp) and 3D DataArray (zeta), no problem:

ds.zeta * ds.temp

<xarray.DataArray (ocean_time: 1325, eta_rho: 602, xi_rho: 677, s_rho: 30)>
dask.array<mul, shape=(1325, 602, 677, 30), dtype=float32, chunksize=(1, 602, 677, 1)>

This returns an object instantly, and the result is lazily evaluated. However, if I just try to add temp to itself,

ds.temp + ds.temp

this fails (eventually) as my medium sized computer runs out of memory, since it starts to evaluate the numbers as if I did a compute() or asked for the values. Note 2*ds.temp or ds.temp**2 is lazily evaluated, and returns an object instantly. Chunk size does not seem to be an issue, as I have tried a number of reasonable choices without success.

Why can't such simple math operations between two large arrays also be lazily evaluated?

This is similar to issue #1503, but specifically deals with requiring the correct order of dimensions for plotting 2d coordinates. An example of the issue is:

```

create an example dataset

dates = pd.date_range('2000-01-01', '2001-12-31', name='time') times = dates - dates[0]

x = np.linspace(0, 10, 101)

h = np.linspace(3, 7, 101) s = np.linspace(0, 1, 51) z = s[:, np.newaxis] * h[np.newaxis, :]

data = (np.sin(x) * np.cos(z)) * np.cos(np.asarray(times.days[:, np.newaxis, np.newaxis]))

create an xarray dataset

ds = xr.Dataset({'data': (('time', 's', 'x'), data)}, {'time':dates, 'x':x, 's':s, 'z':(('s', 'x'), z)}) ds.coords['zt'] = ds.z.transpose() `` plotting withzandzt` should give the same results, since the dimensions are clearly labeled and other than the order of the dimensions, they are exactly the same.

This works: ds.data[0].plot(x='x', y='z')

This does not: ds.data[0].plot(x='x', y='zt')

Broadcasting, e.g., does not give a similar error. Specifically, this is true: np.allclose((ds.data * ds.zt).values, (ds.data * ds.z).values)

• [x] Fully documented, including whats-new.rst for all changes and api.rst for new API

This is an example of reading and visualizing ROMS model output. A highlight is adding a lazily calculated vertical coordinate. This notebook will need an example ROMS output file to run, called ROMS_example.nc. This will be uploaded to the data repository.

I revised the ROMS example notebook to use xr.tutorial.open_dataset, and made some additions to the documentation in the notebook for clarity.

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

Create a new container in variables for CF convention variables. E.g.,

ds.salt.cf.x ds.salt.cf.lon # (maybe also?)

which would pass through ds.salt.x_rho (the x coordinate), or ds.salt.lon_rho.

Create computed coordinates in a lazy way, given a user defined equation for defining these coordinates.

E.g.,

my_depths(params): # do the calculation return depths

Then, link ds.salt.z to my_depths, such that

ds.salt.z

returns the depths, calculated on request. Or better, an object that returns values only upon indexing.

See, e.g., The equations: https://github.com/hetland/octant/blob/master/octant/depths.py

One way to do the lazy calculations just in python https://github.com/hetland/octant/blob/68f51005480c8cb9f83098b0f751b4e8086c11b5/octant/roms.py#L244

This is mainly for plotting in MPL, which does not like NaNs.