It is pretty simply to do with xarray.apply_ufunc.
import gsw
import xarray as xr
import numpy as np
import dask.array as dsa
# define some input data
shape = (100, 1000)
chunks = (100, 200)
sp = xr.DataArray(dsa.full(shape, 35., chunks=chunks), dims=['time', 'depth'])
p = xr.DataArray(np.arange(shape[1]), dims=['depth'])
lon = 0
lat = 45
# lazily apply gsw function
sa = xr.apply_ufunc(gsw.SA_from_SP, sp, p, lon, lat,
dask='parallelized', output_dtypes=[sp.dtype])
# compute in parallel with dask
sa.compute()
A further step to make this even easier would be to have gsw implement NEP 18 and accept dask or xarray arrays as inputs. This would require some changes to gsw itself.