issues

What happened:

A file that isn't obviously malformed or corrupted that I try to load using opendap generates an esoteric IndexError from dask.

What you expected to happen:

For the file to load successfully.

Minimal Complete Verifiable Example:

python url = "http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.HINDCAST/.MONTHLY/.tref/dods" ds_test = xr.open_dataset(url, decode_times=False) print(ds_test) ds_test.load()

This yields the following:

```python <xarray.Dataset> Dimensions: (L: 10, M: 24, S: 348, X: 360, Y: 181) Coordinates: * M (M) float32 1.0 2.0 3.0 4.0 5.0 6.0 ... 20.0 21.0 22.0 23.0 24.0 * X (X) float32 0.0 1.0 2.0 3.0 4.0 ... 355.0 356.0 357.0 358.0 359.0 * L (L) float32 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 * S (S) float32 264.0 265.0 266.0 267.0 ... 608.0 609.0 610.0 611.0 * Y (Y) float32 -90.0 -89.0 -88.0 -87.0 -86.0 ... 87.0 88.0 89.0 90.0 Data variables: tref (S, L, M, Y, X) float32 ... Attributes: Conventions: IRIDL

IndexError Traceback (most recent call last) ~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/backends/netCDF4_.py in _getitem(self, key) 84 original_array = self.get_array(needs_lock=False) ---> 85 array = getitem(original_array, key) 86 except IndexError:

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/backends/common.py in robust_getitem(array, key, catch, max_retries, initial_delay) 53 try: ---> 54 return array[key] 55 except catch:

netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Variable.getitem()

netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Variable._get()

IndexError: index exceeds dimension bounds

During handling of the above exception, another exception occurred:

IndexError Traceback (most recent call last) <ipython-input-9-23bd037d898e> in <module> 2 ds_test = xr.open_dataset(url, decode_times=False) 3 print(ds_test) ----> 4 ds_test.load()

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/core/dataset.py in load(self, **kwargs) 664 for k, v in self.variables.items(): 665 if k not in lazy_data: --> 666 v.load() 667 668 return self

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/core/variable.py in load(self, kwargs) 379 self._data = as_compatible_data(self._data.compute(kwargs)) 380 elif not hasattr(self._data, "array_function"): --> 381 self._data = np.asarray(self._data) 382 return self 383

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/numpy/core/_asarray.py in asarray(a, dtype, order) 81 82 """ ---> 83 return array(a, dtype, copy=False, order=order) 84 85

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/core/indexing.py in array(self, dtype) 675 676 def array(self, dtype=None): --> 677 self._ensure_cached() 678 return np.asarray(self.array, dtype=dtype) 679

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/core/indexing.py in _ensure_cached(self) 672 def _ensure_cached(self): 673 if not isinstance(self.array, NumpyIndexingAdapter): --> 674 self.array = NumpyIndexingAdapter(np.asarray(self.array)) 675 676 def array(self, dtype=None):

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/numpy/core/_asarray.py in asarray(a, dtype, order) 81 82 """ ---> 83 return array(a, dtype, copy=False, order=order) 84 85

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/core/indexing.py in array(self, dtype) 651 652 def array(self, dtype=None): --> 653 return np.asarray(self.array, dtype=dtype) 654 655 def getitem(self, key):

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/numpy/core/_asarray.py in asarray(a, dtype, order) 81 82 """ ---> 83 return array(a, dtype, copy=False, order=order) 84 85

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/core/indexing.py in array(self, dtype) 555 def array(self, dtype=None): 556 array = as_indexable(self.array) --> 557 return np.asarray(array[self.key], dtype=None) 558 559 def transpose(self, order):

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/backends/netCDF4_.py in getitem(self, key) 71 def getitem(self, key): 72 return indexing.explicit_indexing_adapter( ---> 73 key, self.shape, indexing.IndexingSupport.OUTER, self._getitem 74 ) 75

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/core/indexing.py in explicit_indexing_adapter(key, shape, indexing_support, raw_indexing_method) 835 """ 836 raw_key, numpy_indices = decompose_indexer(key, shape, indexing_support) --> 837 result = raw_indexing_method(raw_key.tuple) 838 if numpy_indices.tuple: 839 # index the loaded np.ndarray

~/miniconda3/envs/ensinoo/lib/python3.7/site-packages/xarray/backends/netCDF4_.py in _getitem(self, key) 93 "your data into memory first by calling .load()." 94 ) ---> 95 raise IndexError(msg) 96 return array 97

IndexError: The indexing operation you are attempting to perform is not valid on netCDF4.Variable object. Try loading your data into memory first by calling .load(). ```

Anything else we need to know?:

There's something specific to this file for sure that's causing this, as I'm able to load (for expediency, a subset of) the file from the xarray docs section on opendap successfully:

```python url = "http://iridl.ldeo.columbia.edu/SOURCES/.OSU/.PRISM/.monthly/.tdmean/[X+]average/dods" ds_test = xr.open_dataset(url, decode_times=False) print(ds_test)

<xarray.Dataset> Dimensions: (T: 1420, Y: 621) Coordinates: * Y (Y) float32 49.916668 49.875 49.833336 ... 24.125 24.083334 * T (T) float32 -779.5 -778.5 -777.5 -776.5 ... 636.5 637.5 638.5 639.5 Data variables: tdmean (T, Y) float64 ... Attributes: Conventions: IRIDL ```

And then calling ds_test.load() works just fine.

Environment:

• [x] Fixes #3287
• [x] Tests added
• [x] Passes isort -rc . && black . && mypy . && flake8
• [x] Fully documented, including whats-new.rst for all changes and api.rst for new API

I've added a check within groupby.__init__ for monotonicity of the coords, which makes my test pass, but it causes three existing tests to fail. Getting the logic right that properly distinguishes among these different cases is escaping me for now, but it does seem like this call to unique_value_groups within groupby.__init__ is where the re-sorting is occurring. Feedback welcome. TIA!

A heads up: @spencerkclark unearthed problems with libnetcdf 4.5.0 that will cause xarray test failures. See https://github.com/Unidata/netcdf4-python/issues/742 and link therein for more info.

Now that CFTimeIndex has been implemented (#1252), one thing that remains to implement is resampling. @shoyer provided a sketch of how to implement it: https://github.com/pydata/xarray/pull/1252#issuecomment-380593243. In the interim, @spencerkclark provided a sketch of a workaround for some use-cases using groupby: https://github.com/pydata/xarray/issues/1270#issuecomment-390973986.

I thought it would be useful to have a new Issue specifically on this topic from which future conversation can continue. @shoyer, does that sketch you provided still seem like a good starting point?

DataArray.mean applied along a particular dimension causes coordinates that are defined in that dimension to be dropped:

```python In [39]: x, y = range(2), range(3) In [40]: arr = xr.DataArray(np.random.random((2,3)), dims=['x', 'y'], coords=dict(x=x, y=y)) In [41]: coord = xr.DataArray(np.random.random((2,3)), dims=['x', 'y'], coords=dict(x=x, y=y)) In [42]: arr = arr.assign_coords(z=coord) In [43]: arr Out[43]: <xarray.DataArray (x: 2, y: 3)> array([[ 0.132368, 0.746242, 0.48783 ], [ 0.12747 , 0.751283, 0.033713]]) Coordinates: * y (y) int64 0 1 2 * x (x) int64 0 1 z (x, y) float64 0.993 0.1031 0.1808 0.2769 0.7237 0.2891

In [44]: arr.mean('x') Out[44]: <xarray.DataArray (y: 3)> array([ 0.129919, 0.748763, 0.260772]) Coordinates: * y (y) int64 0 1 2 ```

We have a use case where we'd like to preserve the coordinates. @spencerkclark came up with the following workaround, which entails converting to a dataset, promoting the coord to a variable, performing the mean, and demoting the original coord back from a variable to a coord:

python def coord_preserving_mean(arr, *args, **kwargs): name = arr.name ds = arr.reset_coords() names = set(arr.coords) - set(arr.dims) ds = ds.mean(*args, **kwargs) return ds.set_coords(names)[name]

This works fine, but I just feel like maybe there's an easier way to do this that we're missing. Any ideas? Thanks in advance.

xref https://github.com/spencerahill/aospy/issues/194

Re: #789. The consensus is that upstream fixes in Pandas are not coming anytime soon, and there is an acute need amongst many xarray users for a workaround in the meantime. There are two separate issues: (1) date-range limitations due to nanosecond precision, and (2) support for non-standard calendars.

@shoyer, @jhamman , @spencerkclark, @darothen, and I briefly discussed offline a potential workaround that I am (poorly) summarizing here, with hope that others will correct/extend my snippet.

The idea is to extend either PeriodIndex or (more involved but potentially more robust) Int64Index, either through subclassing or composition, to implement all of the desired functionality: slicing, resampling, groupby, and serialization.

For reference, @spencerkclark nicely summarized the limitations of PeriodIndex and the netCDF4.datetime objects, which are often used as workarounds currently: https://github.com/spencerahill/aospy/issues/98#issuecomment-256043833

Just want to clarify if below is the intended behavior re: selecting on dims w/ vs. w/out coords. I'm on v0.9.1:

In short, selecting a single value of a dimension causes the dimension to be dropped only if it doesn't have a coordinate. Is that intentional? It's counter-intuitive to me but not a huge deal. Simple example:

```python In [156]: ds=xr.DataArray([1,2], name='test').to_dataset().assign_coords(bounds=[0,1])

In [161]: ds Out[161]: <xarray.Dataset> Dimensions: (bounds: 2, dim_0: 2) Coordinates: * bounds (bounds) int64 0 1 Dimensions without coordinates: dim_0 Data variables: test (dim_0) int64 1 2

In [163]: ds.isel(bounds=0).drop('bounds') Out[163]: <xarray.Dataset> Dimensions: (dim_0: 2) Dimensions without coordinates: dim_0 Data variables: test (dim_0) int64 1 2

In [164]: ds.isel(dim_0=0).drop('dim_0')

ValueError Traceback (most recent call last) <ipython-input-164-902de222762c> in <module>() ----> 1 ds.isel(dim_0=0).drop('dim_0')

/Users/shill/Dropbox/miniconda3/lib/python3.5/site-packages/xarray/core/dataset.py in drop(self, labels, dim) 1888 labels = [labels] 1889 if dim is None: -> 1890 return self._drop_vars(labels) 1891 else: 1892 try:

/Users/shill/Dropbox/miniconda3/lib/python3.5/site-packages/xarray/core/dataset.py in _drop_vars(self, names) 1899 1900 def _drop_vars(self, names): -> 1901 self._assert_all_in_dataset(names) 1902 drop = set(names) 1903 variables = OrderedDict((k, v) for k, v in iteritems(self._variables)

/Users/shill/Dropbox/miniconda3/lib/python3.5/site-packages/xarray/core/dataset.py in _assert_all_in_dataset(self, names, virtual_okay) 1867 bad_names -= self.virtual_variables 1868 if bad_names: -> 1869 raise ValueError('One or more of the specified variables ' 1870 'cannot be found in this dataset') 1871

ValueError: One or more of the specified variables cannot be found in this dataset ``` cc @spencerkclark xref https://github.com/spencerahill/aospy/issues/137#issuecomment-277326319

I.e. there is no from . import ufuncs in __init__.py. Is this deliberate?

If it was included, one could access them after import xarray as xr via e.g. xr.ufuncs.sin rather than having to do separately e.g. from xarray.ufuncs import sin.

Sorry if I'm missing something obvious here.

• [x] closes #1216
• [ ] tests added / passed
• [ ] passes git diff upstream/master | flake8 --diff
• [ ] whatsnew entry

Let me know if you want tests and/or what's new; wasn't sure since this is so small.

Similar to numpy's numpy.testing module, I would find it useful for xarray to expose some of its testing methods, particularly those in its base TestCase.

I imagine other folks whose packages and/or modules derive heavily from xarray would find this useful too. For now I am just copy-pasting xarray's code for these methods into my testing modules.

I think this is another unintended consequence of #648. Consider the following case:

python In [20]: arr = xr.DataArray(range(3), dims=['abc']) In [21]: new_coord = xr.DataArray([1,2,3], dims=['abc'], coords=[[1,2,3]]) In [22]: arr['abc'] = new_coord In [23]: arr Out[24]: <xarray.DataArray (abc: 3)> array([0, 1, 2]) Coordinates: * abc (abc) float64 nan 1.0 2.0

Note the nan. Before #648, this worked, in that arr's coordinates would consist of new_coord.

The use case: we have some data defined on the edges of pressure levels in an atmospheric model, and other data defined at the center of the pressure levels. In order to perform calculations involving both kinds of data, we average the edge-defined data (i.e. 0.5*(value at top edge + value at bottom edge)) to get the value at the level centers. But the resulting DataArray still has as its coord (from xarray's perspective, that is) the level edges, and so we replace that coord with the DataArray of the level centers.

A workaround would be arr['abc'].values = new_coord.values, but then any other metadata associated with the original coord is retained, which is not good.

Somewhat involved and not sure I described clearly, so let me know if clarification needed. Also I vaguely suspect there's a cleaner way of doing this in the first place. Thanks!

Closes #724

Before #648, passing Dataset.rename a name_dict whose keys and values were identical did not raise an exception. Now, however, it raises a ValueError:

``` python arr = xr.DataArray(range(2), name='arrname') ds = arr.to_dataset() ds.rename({'dim_0':'dim_0'})

ValueError Traceback (most recent call last) <ipython-input-15-a5b851f9fb10> in <module>() ----> 1 ds.rename({'dim_0':'dim_0'})

/Users/spencerahill/anaconda/lib/python2.7/site-packages/xarray/core/dataset.pyc in rename(self, name_dict, inplace) 1245 "variable in this dataset" % k) 1246 if v in self: -> 1247 raise ValueError('the new name %r already exists' % v) 1248 1249 variables = OrderedDict()

ValueError: the new name 'dim_0' already exists ```

This is easy enough to handle with a try/except clause in my own code, but it would be nice (for me at least) for rename to not raise anything for these cases. Since the result is simply the original Dataset unchanged, this could be implemented by replacing line 1246 with if v in self and k != v:

The use case is that we have data coming from multiple sources, often with differing internal names for coordinates and variables, and we're programmatically forcing them to have consistent names via rename. Sometimes, a piece of data already has the name that we are after, hence the non-rename rename. Thanks!

DataArray.diff, when applied to a coordinate, just returns the coordinate itself with one index truncated. Consider the following:

``` In [5]: arr = xray.DataArray(range(0, 20, 2), dims=['lon'], coords=[range(10)]) In [6]: arr.diff('lon') <xray.DataArray (lon: 9)> array([2, 2, 2, 2, 2, 2, 2, 2, 2]) Coordinates: * lon (lon) int64 1 2 3 4 5 6 7 8 9

In [7]: arr['lon'].diff('lon') <xray.DataArray 'lon' (lon: 9)> array([1, 2, 3, 4, 5, 6, 7, 8, 9]) Coordinates: * lon (lon) int64 1 2 3 4 5 6 7 8 9 ```

Is this the intended behavior? The documentation doesn't mention anything about this, and its counter-intuitive, so I'm wondering if its a bug instead. Even if it is intended, I personally would like to be able to use a coordinate array's diff on itself to get its spacing, e.g. for use as the denominator in finite differencing approximations to derivatives.

Thanks!

Most of the netCDF data I work with stores, in addition to the coordinates themselves, the bounds of each coordinate value. Often these bounds are stored as arrays with shape Nx2, where N is the number of points for that coordinate. For example:

``` $ ncdump -c /archive/Spencer.Hill/am3/am3clim_hurrell/gfdl.ncrc2-intel-prod-openmp/pp/atmos/ts/monthly/1yr/atmos.201001-201012.t_surf.nc netcdf atmos.201001-201012.t_surf { dimensions: time = UNLIMITED ; // (12 currently) lat = 90 ; bnds = 2 ; lon = 144 ; variables: double average_DT(time) ; average_DT:long_name = "Length of average period" ; average_DT:units = "days" ; average_DT:missing_value = 1.e+20 ; average_DT:_FillValue = 1.e+20 ; double average_T1(time) ; average_T1:long_name = "Start time for average period" ; average_T1:units = "days since 1980-01-01 00:00:00" ; average_T1:missing_value = 1.e+20 ; average_T1:_FillValue = 1.e+20 ; double average_T2(time) ; average_T2:long_name = "End time for average period" ; average_T2:units = "days since 1980-01-01 00:00:00" ; average_T2:missing_value = 1.e+20 ; average_T2:_FillValue = 1.e+20 ; double lat(lat) ; lat:long_name = "latitude" ; lat:units = "degrees_N" ; lat:cartesian_axis = "Y" ; lat:bounds = "lat_bnds" ; double lat_bnds(lat, bnds) ; lat_bnds:long_name = "latitude bounds" ; lat_bnds:units = "degrees_N" ; lat_bnds:cartesian_axis = "Y" ; double lon(lon) ; lon:long_name = "longitude" ; lon:units = "degrees_E" ; lon:cartesian_axis = "X" ; lon:bounds = "lon_bnds" ; double lon_bnds(lon, bnds) ; lon_bnds:long_name = "longitude bounds" ; lon_bnds:units = "degrees_E" ; lon_bnds:cartesian_axis = "X" ; float t_surf(time, lat, lon) ; t_surf:long_name = "surface temperature" ; t_surf:units = "deg_k" ; t_surf:valid_range = 100.f, 400.f ; t_surf:missing_value = 1.e+20f ; t_surf:_FillValue = 1.e+20f ; t_surf:cell_methods = "time: mean" ; t_surf:time_avg_info = "average_T1,average_T2,average_DT" ; t_surf:interp_method = "conserve_order2" ; double time(time) ; time:long_name = "time" ; time:units = "days since 1980-01-01 00:00:00" ; time:cartesian_axis = "T" ; time:calendar_type = "JULIAN" ; time:calendar = "JULIAN" ; time:bounds = "time_bounds" ; double time_bounds(time, bnds) ; time_bounds:long_name = "time axis boundaries" ; time_bounds:units = "days" ; time_bounds:missing_value = 1.e+20 ; time_bounds:_FillValue = 1.e+20 ;

// global attributes: :filename = "atmos.201001-201012.t_surf.nc" ; :title = "am3clim_hurrell" ; :grid_type = "mosaic" ; :grid_tile = "1" ; :comment = "pressure level interpolator, version 3.0, precision=double" ; :history = "fregrid --input_mosaic atmos_mosaic.nc --input_file 20100101.atmos_month --interp_method conserve_order2 --remap_file .fregrid_remap_file_144_by_90 --nlon 144 --nlat 90 --scalar_field (please see the field list in this file)" ; :code_version = "$Name: fre-nctools-bronx-7 $" ; data:

lat = -89, -87, -85, -83, -81, -79, -77, -75, -73, -71, -69, -67, -65, -63, -61, -59, -57, -55, -53, -51, -49, -47, -45, -43, -41, -39, -37, -35, -33, -31, -29, -27, -25, -23, -21, -19, -17, -15, -13, -11, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, 77, 79, 81, 83, 85, 87, 89 ;

lon = 1.25, 3.75, 6.25, 8.75, 11.25, 13.75, 16.25, 18.75, 21.25, 23.75, 26.25, 28.75, 31.25, 33.75, 36.25, 38.75, 41.25, 43.75, 46.25, 48.75, 51.25, 53.75, 56.25, 58.75, 61.25, 63.75, 66.25, 68.75, 71.25, 73.75, 76.25, 78.75, 81.25, 83.75, 86.25, 88.75, 91.25, 93.75, 96.25, 98.75, 101.25, 103.75, 106.25, 108.75, 111.25, 113.75, 116.25, 118.75, 121.25, 123.75, 126.25, 128.75, 131.25, 133.75, 136.25, 138.75, 141.25, 143.75, 146.25, 148.75, 151.25, 153.75, 156.25, 158.75, 161.25, 163.75, 166.25, 168.75, 171.25, 173.75, 176.25, 178.75, 181.25, 183.75, 186.25, 188.75, 191.25, 193.75, 196.25, 198.75, 201.25, 203.75, 206.25, 208.75, 211.25, 213.75, 216.25, 218.75, 221.25, 223.75, 226.25, 228.75, 231.25, 233.75, 236.25, 238.75, 241.25, 243.75, 246.25, 248.75, 251.25, 253.75, 256.25, 258.75, 261.25, 263.75, 266.25, 268.75, 271.25, 273.75, 276.25, 278.75, 281.25, 283.75, 286.25, 288.75, 291.25, 293.75, 296.25, 298.75, 301.25, 303.75, 306.25, 308.75, 311.25, 313.75, 316.25, 318.75, 321.25, 323.75, 326.25, 328.75, 331.25, 333.75, 336.25, 338.75, 341.25, 343.75, 346.25, 348.75, 351.25, 353.75, 356.25, 358.75 ;

time = 10973.5, 11003, 11032.5, 11063, 11093.5, 11124, 11154.5, 11185.5, 11216, 11246.5, 11277, 11307.5 ; } ```

These 2-D bounding arrays lead to the "Buffer has wrong number of dimensions" error in #665. In the case of #665, only the time coordinate has this 2-D bounds array; here other coordinates (namely lat and lon) have it as well.

Conceptually, these bound arrays represent coordinates, but when read in as a Dataset, they become variables, not coordinates. Perhaps this is part of the problem?