Hello,

When I run this line:

ds = cat["en4"].to_dask()

sometimes, no problem; however, sometimes, there's an error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-24-dbd9e41d9ee1> in <module>()
----> 1 ds = cat["en4"].to_dask()
      2 print("Size of the dataset:", ds.nbytes/1e9,"Gb")
      3 ds

7 frames
/usr/local/lib/python3.7/dist-packages/xarray/backends/plugins.py in get_backend(engine)
    133                 f"backends {installed_engines}. Consider explicitly selecting one of the "
    134                 "installed engines via the ``engine`` parameter, or installing "
--> 135                 "additional IO dependencies, see:\n"
    136                 "http://xarray.pydata.org/en/stable/getting-started-guide/installing.html\n"
    137                 "http://xarray.pydata.org/en/stable/user-guide/io.html"

ValueError: unrecognized engine zarr must be one of: ['netcdf4', 'scipy', 'store']

I tried !pip install xarray[complete] as https://github.com/pydata/xarray/issues/5395 suggested, but it still didn't work.

Hello!

We have three doubts:

1. To compute Conservative Temperature, should we use gsw.CT_from_t(ds.SA, ds.temperature, ds.P) , giving the pressure at different levels or gsw.CT_from_pt(ds.SA, ds.temperature), where reference pressure is taken at 0 dbar?
2. To apply gsw.geo_strf_dyn_height which pressure reference should we use? Also I would like to know which dimensions are used here. In the docs TEOS-10 it first defines pref as 10 dbar but later in the example it applies 1000 without any units.
3. How should we pass from steric height anomaly to steric height?

These doubts are better specified in the code found at the end.

We think we are not really understanding how to use the pressure reference coherently as the value of the steric height anomaly we are obtaining is of the order of metres and we are expecting values in the order of millimetres.

# Necessary imports
import numpy as np
import xarray as xr
import dask
import matplotlib.pyplot as plt
import gsw
from intake import open_catalog

# Open dataset:
catalog_url = 'https://raw.githubusercontent.com/obidam/ds2-2022/main/ds2_data_catalog.yml'
cat = open_catalog(catalog_url)
ds = cat["en4"].to_dask()
# Select only first 34 depth levels
ds = ds.sel(depth=slice(ds.depth[0],ds.depth[34]))

# Treat data:
## Pass depth values to negative (positive direction upwards).
ds['depth'] = - ds['depth']
## Pass temperature values to degrees celsius.
ds['temperature'] = ds['temperature'] - 273.15
## Calculate pressure values from depth + latitude.
### Then make it coincide with Dataset dimensions so it can be 
### applied to other gsw functions:
ds['P'] = gsw.p_from_z(ds.depth, ds.lat).transpose('lat', 'depth') .expand_dims({'lon' : ds.lon}).where(
                            ds['temperature'] == ds['temperature']).transpose('time', 'depth', 'lat', 'lon') 

## Calculate Absolute Salinity and Conservative Temperature.
ds['SA'] = gsw.SA_from_SP(ds.salinity,ds.P,ds.lon,ds.lat)
##!!! Here we have a doubt among which one to use for CT:
ds['CT'] = gsw.CT_from_t(ds.SA, ds.temperature, ds.P)
## --> or should we use this one?:
##ds['CT'] = gsw.CT_from_pt(ds.SA, ds.temperature)


## Average among years. Put pressure dimension axis as first axis for dynamic height.
ds = ds.groupby('time.year').mean('time').transpose('depth', 'year', 'lat', 'lon','bnds')

## Apply the dynamic height anomaly parallelized function and divide by the gravity acceleration to obtain steric height anomaly.
## We also have a doubt if the reference  should  be equal to 1000 ¿?¿?
g_0 = 9.7963 # m/s^2 
ds['st_height_anom'] = xr.apply_ufunc( gsw.geo_strf_dyn_height, ds.SA, ds.CT, ds.P, 1000, dask='parallelized', output_dtypes=[ds.SA.dtype])/g_0

# Testing for 4 years, should take less than 1 min.
ds = ds.sel(year=slice(2015,2019))
(ds['st_height_anom'].isel(depth=0).mean(dim=('lat','lon'))).plot()