issue_comments

I think you're right, this is nearly the same issue as #2850. Your dataset has an attribute named _NCProperties which is one of the hidden/forbidden attributes.

If I remove the attribute: del dsSel.attrs['_NCProperties']

dsSel.to_netcdf(outnc) then works for me with no errors.

You can also delete the attribute on ds instead and subsequent selections should also work.

As an aside, if I run ncdump -h https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/dodsC/birdhouse/cccs_portal/indices/Final/BCCAQv2/tg_mean/YS/rcp26/simulations/BCCAQv2+ANUSPLIN300_bcc-csm1-1_historical+rcp26_r1i1p1_1950-2100_tg_mean_YS.nc "_NCProperties" doesn't show up in the attributes, but I assume this is because ncdump hides it and not because it gets added by xarray or netcdf4. If that's not the case it could be a bug in xarray or netcdf4.

Sounds great. Thanks for making this a smooth and well documented process for new contributors

I have to say as someone who is probably an average user, inconsistencies between related projects, like xarray/pandas or matplotlib/seaborn, drive me nuts. But I also agree that missing= is more descriptive, and if we were starting from scratch I would totally support that. So I will defer to the maintainers here.

I can reproduce OP when the coordinate is a datetime64, but the problem does not occur with int64s (or float64s).

datetime coordinate: ``` In [1]: import xarray

In [2]: xarray.version Out[2]: '0.11.2'

In [3]: import datetime as dt

In [4]: da = xarray.DataArray([0], coords=[('time', [dt.datetime(1979, 1, 1)])])

In [5]: da.isel(time=0).time.count() Out[5]: <xarray.DataArray 'time' ()> array(-1) Coordinates: time datetime64[ns] 1979-01-01 ```

int coordinate:

``` In [6]: da = xarray.DataArray([0], coords=[('time', [0])])

In [7]: da.isel(time=0).time.count() Out[7]: <xarray.DataArray 'time' ()> array(1) Coordinates: time int64 0 ```

I think this might be a problem with bottleneck? My interpretation of _create_nan_agg_method in xarray/core/ops.py is that it may use bottleneck to get the mean unless you pass skipna=False or specify multiple axes. And,

```python In [2]: import bottleneck In [3]: bottleneck.version Out[3]: '1.2.0'

In [6]: bottleneck.nanmean(ds.var167.data) Out[6]: 261.6441345214844 ```

Forgive me if I'm wrong, I'm still a bit new.

Also on macOS, and I can reproduce.

Using python 2.7.11, xarray 0.9.1, dask 0.14.1 installed through Anaconda. I get the same results with xarray 0.9.1-38-gc0178b7 from GitHub.

```python In [3]: ds = xarray.open_dataset('ERAIN-t2m-1983-2012.seasmean.nc')

In [4]: ds.var167.mean() Out[4]: <xarray.DataArray 'var167' ()> array(261.6441345214844, dtype=float32) ```

Curiously, I get the right results with skipna=False...

python In [10]: ds.var167.mean(skipna=False) Out[10]: <xarray.DataArray 'var167' ()> array(278.6246643066406, dtype=float32)

... or by specifying coordinates to average over:

python In [5]: ds.var167.mean(('time', 'lat', 'lon')) Out[5]: <xarray.DataArray 'var167' ()> array(278.6246643066406, dtype=float32)