Non-zeta-like parameters fail to .at() with v0.0.9 Release

[acrosby]

I am getting the following error trying to access non-zeta (and non-depth avg velocity) parameters from SGrid datasets. This example by @rsignell-usgs (https://github.com/NOAA-ORR-ERD/gridded/blob/master/examples/trajectory_test.ipynb) indicates that the functionality should work, but I have tested with the example BORA dataset and get the same errors.

/usr/local/lib/python2.7/dist-packages/gridded/variable.pyc in at(self, points, time, units, extrapolate, _hash, _mem, _auto_align, **kwargs)
    424         order = self.dimension_ordering
    425         if order[0] == 'time':
--> 426             value = self._time_interp(pts, time, extrapolate, _mem=_mem, _hash=_hash, **kwargs)
    427         elif order[0] == 'depth':
    428             value = self._depth_interp(pts, time, extrapolate, _mem=_mem, _hash=_hash, **kwargs)

/usr/local/lib/python2.7/dist-packages/gridded/variable.pyc in _time_interp(self, points, time, extrapolate, slices, **kwargs)
    476         if time == self.time.min_time or (extrapolate and time < self.time.min_time):
    477             # min or before
--> 478             return val_func(points, time, extrapolate, slices=(0,), ** kwargs)
    479         elif time == self.time.max_time or (extrapolate and time > self.time.max_time):
    480             return val_func(points, time, extrapolate, slices=(-1,), **kwargs)

/usr/local/lib/python2.7/dist-packages/gridded/variable.pyc in _depth_interp(self, points, time, extrapolate, slices, **kwargs)
    508         else:
    509             val_func = self._time_interp
--> 510         indices, alphas = self.depth.interpolation_alphas(points, time, self.data.shape[1:], kwargs.get('_hash', None))
    511         if indices is None and alphas is None:
    512             # all particles are on surface

/usr/local/lib/python2.7/dist-packages/gridded/depth.pyc in interpolation_alphas(self, points, time, data_shape, _hash)
    342         index and depth_index+1. If both values are None, then all particles are on the surface layer.
    343         '''
--> 344         underwater = points[:, 2] > -self.zeta.at(points, time)
    345         if len(np.where(underwater)[0]) == 0:
    346             return None, None

IndexError: index 2 is out of bounds for axis 1 with size 2

I have a similar problem with UGrid datasets, but different error. Both issues seem to be occurring in the part of the code that is trying to deal with depth. Which explains the parameters that are being impacted.

/usr/local/lib/python2.7/dist-packages/gridded/variable.pyc in at(self, points, time, units, extrapolate, _hash, _mem, _auto_align, **kwargs)
    424         order = self.dimension_ordering
    425         if order[0] == 'time':
--> 426             value = self._time_interp(pts, time, extrapolate, _mem=_mem, _hash=_hash, **kwargs)
    427         elif order[0] == 'depth':
    428             value = self._depth_interp(pts, time, extrapolate, _mem=_mem, _hash=_hash, **kwargs)

/usr/local/lib/python2.7/dist-packages/gridded/variable.pyc in _time_interp(self, points, time, extrapolate, slices, **kwargs)
    476         if time == self.time.min_time or (extrapolate and time < self.time.min_time):
    477             # min or before
--> 478             return val_func(points, time, extrapolate, slices=(0,), ** kwargs)
    479         elif time == self.time.max_time or (extrapolate and time > self.time.max_time):
    480             return val_func(points, time, extrapolate, slices=(-1,), **kwargs)

/usr/local/lib/python2.7/dist-packages/gridded/variable.pyc in _depth_interp(self, points, time, extrapolate, slices, **kwargs)
    508         else:
    509             val_func = self._time_interp
--> 510         indices, alphas = self.depth.interpolation_alphas(points, time, self.data.shape[1:], kwargs.get('_hash', None))
    511         if indices is None and alphas is None:
    512             # all particles are on surface

TypeError: interpolation_alphas() takes at most 4 arguments (5 given)

Are there canonical SGrid and UGrid datasets that gridded is tested against that I can use to rule out issues with the datasets grid/topology/depth identification?

[ChrisBarker-NOAA]

I'll let @jay-hennen try to debug this, but as for the canonical datasets, no. But we do have the ones we happen to have tested with. @jay-hennen should be able to point you to those.

But it's a good idea -- we should have a good set for tests and demonstrations. I"d like to see both:

• a full set of tiny examples in with the test code.
• a full set of "real" scale examples -- those shouldn't be in the repo, but we can host them on a public-facing server. Or maybe there's a public server somewhere we could host them. Or openDAP urls -- but that feels a bit too ephemeral.

We'd love anyone to help find a comprehensive set to test with.

[acrosby]

@ChrisBarker-NOAA I'm happy if to host some data in our AWS S3 buckets with public http access or at least until a better solution can be found.

[jay-hennen]

Can you provide more metadata about the dataset or variable you're trying to open (shape? level or zeta depth?) Also, the function call you used to create the Variable in the first place

What's happening in the first case seems to be a misinterpretation of the data. It's using a zeta-style depth interpolator despite the data being level-depth (right?)

The second might be a straight bug, but I'd need to check the metadata to be sure...

[acrosby]

For the first SGrid related error (also same issue using the WHOI BORA dataset found in the example referenced above), Metadata for dataset: https://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/TBOFS/MODELS/201801/nos.tbofs.fields.f006.20180116.t00z.nc.html

f = "http://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/TBOFS/MODELS/201801/nos.tbofs.fields.f006.20180116.t00z.nc"

nc = gridded.Dataset(f)
bbox = get_bbox(nc)

var = "salt" 
#var = "u"
#var = "zeta"  # works

p = nc.variables[var] 
time = p.time.get_time_array()
rlon, rlat = np.meshgrid(np.linspace(bbox[0], bbox[2], 1000),
                         np.linspace(bbox[1], bbox[3], 1000))
rp = p.at(np.vstack((rlon.ravel(),
                         rlat.ravel())).T, 
              time[0]).reshape(rlon.shape)

For the second error message (UGrid, possibly related to #19, but in general the api goes seem to know when vars are nodal or not), Metadata: http://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/NWGOFS/MODELS/201801/nos.nwgofs.fields.f006.20180116.t03z.nc.html

f = "http://opendap.co-ops.nos.noaa.gov/thredds/dodsC/NOAA/NWGOFS/MODELS/201801/nos.nwgofs.fields.f006.20180116.t03z.nc"

nc = gridded.Dataset(f)
bbox = get_bbox(nc)

var = "salinity" 
#var = "u"
#var = "zeta"  # works
#var = "ua"  # works

p = nc.variables[var] 
time = p.time.get_time_array()
rlon, rlat = np.meshgrid(np.linspace(bbox[0], bbox[2], 1000),
                         np.linspace(bbox[1], bbox[3], 1000))
rp = p.at(np.vstack((rlon.ravel(),
                         rlat.ravel())).T, 
              time[0]).reshape(rlon.shape)

[jay-hennen]

Looks like in the SGrid case, you aren't using 3D points in the .at call, if I'm reading the code right.

The UGrid problem is that I think that the depth coordinate system described in the metadata isn't actually supported yet. I'll have to look deeper at the siglay and siglev variables to be sure.

[ChrisBarker-NOAA]

Looks like in the SGrid case, you aren't using 3D points in the .at call, if I'm reading the code right.

if that is the case, let's make sure it gives a better error msg! Or was it giving a 2D result OK, just not what you were expecting?

[jay-hennen]

Yeah, it should be able to give a better message than that. It should only give it in the _depth_interp function, since that's when 3D points become required

[acrosby]

So am I understanding correctly that the points input should have a third coordinate for depth? What does the depth kwarg do? --or perhaps the docstring is not up-to-date?

EDIT: Yes this does appear to work for me, very cool! I guess the docstring and example need to be updated?

[ChrisBarker-NOAA]

PRs accepted :-)

[acrosby]

@jay-hennen Did you find out if this is the case?

The UGrid problem is that I think that the depth coordinate system described in the metadata isn't actually supported yet. I'll have to look deeper at the siglay and siglev variables to be sure.

Actually, with the warnings being reported in jupyter, I can see that this is probably the case:

/usr/local/lib/python2.7/dist-packages/gridded/depth.py:507: RuntimeWarning: Unable to automatically determine depth system so reverting to surface-only mode. Please verify the index of the surface is correct.
  warnings.warn('''Unable to automatically determine depth system so reverting to surface-only mode. Please verify the index of the surface is correct.''', RuntimeWarning)