surface dataset: gradient calculation

[apatlpo]

I'm pretty sure this is a usage mistake, but I am having difficulties computing a gradient:

>> print(ds.eta)
<xarray.DataArray 'eta' (time: 24, Nc: 225368)>
dask.array<getitem, shape=(24, 225368), dtype=float64, chunksize=(24, 2000), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 2013-07-01T00:01:00 ... 2013-07-01T23:00:00
Dimensions without coordinates: Nc
Attributes:
    location:   face
    long_name:  Sea surface elevation
    mesh:       suntans_mesh
    units:      m

>> ds.suntans.calc_grad(ds.eta)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[10], line 1
----> 1 ds.suntans.calc_grad(ds.eta)

File ~/nwa/nwatools/uplot.py:214, in Plot.calc_grad(self, phi, k)
    211 mask = ne == self._FillValue
    212 ne[mask]=0
--> 214 Gn_phi = _GradientAtFace(phi,ne,k)
    215 Gn_phi[mask]=0
    217 Gx_phi = Gn_phi * self.n1[ne] * self.DEF * self.df[ne]

File ~/nwa/nwatools/uplot.py:191, in Plot.calc_grad.<locals>._GradientAtFace(phi, jj, k)
    189 grad1 = self.grad[:,0]
    190 grad2 = self.grad[:,1]
--> 191 nc1 = grad1[jj]
    192 nc2 = grad2[jj]
    194 # check for edges (use logical indexing)

IndexError: index -999999 is out of bounds for axis 0 with size 679391

Any pointers as to how to fix this are welcome

[apatlpo]

For the records, here is the correct way to compute the gradient:

ds = xr.open_zarr(zarr)

# manually negate 999999 values
ds["cells"] = ds.cells.where( ds.cells!=999999, other=-999999 )
# manually add Nk to suntans accessor
ds.suntans.Nk = np.ones(ds.Nc.size)

# compute gradient on a temporal snapshot provided as a numpy array
ds.suntans.calc_grad(ds.eta.isel(time=1))