issue_comments

Yes, seems like a redirect issue. The URL is fine.

Thanks for reminding me, stalte bot, this is still relevant to me.

I think your proposed solution @shoyer would solve this issue perfectly: expose align's join parameter in the concat function - and please also in open_mfdataset. Would you like a PR for this?

The original issue of this thread is that you sometimes might want to disable alignment checks for coordinates other than the concat_dim and only check for same dimensions and dimension shapes.

When you xr.merge with join='exact', it still checks for alignment (see https://github.com/pydata/xarray/pull/1330#issuecomment-302711852), but does not join the coordinates if they are not aligned. This behavior (not joining) is also included in what @rabernat envisioned here, but his suggestion goes beyond that: you don't even load coordinate values from all but the first dataset and just blindly trust that they are aligned.

So xr.open_mfdataset(join='exact', coords='minimal') does not fix this issue here, I think.

I'm confused. It seems that I did exactly this a while ago

@fmaussion don't be. You are perfectly right, we worked on the same line of code, only that I was about 50 PRs earlier, hence the merge conflict.

I'll close this.

That single AppVeyor CI fail is because the py27 32bit conda environment cannot be created: ResolvePackageNotFound: - cftime

mind freshening up this PR

I deleted my fork in the meantime. I opened a new PR at #2096 and fixed the merge conflict there.

Can you add a brief note on the bug fix in whats-new?

Done.

I can see netcdftime is casting the time numbers to float64:

https://github.com/Unidata/netcdftime/blob/e745434547de728d53cde6316cae75b62db102e2/netcdftime/_netcdftime.pyx#L241

like I also suggested for https://github.com/pydata/xarray/pull/1863#issuecomment-361089977.

How about we do that, too?

Seems like this issue only exists on Windows. See https://ci.appveyor.com/project/shoyer/xray/build/1.0.3782/job/yjodrbenev1tw5y7#L4686

(The Travis failure is due to a failing build.)

Can't this be closed?

I found out why the issue is absent in the current implementation. Please see https://github.com/pydata/xarray/pull/1863#issuecomment-361091073. That PR adds the test you asked for @shoyer. It passes in the current version and fails in v0.10.0.

Store _NS_PER_TIME_DELTA values as int, then numpy will do the casting.

Haha, OK, this is actually what you implemented in https://github.com/pydata/xarray/commit/50b0a69a7aa0fb7ac3afb28e7bd971cf08055f99:

https://github.com/pydata/xarray/blob/50b0a69a7aa0fb7ac3afb28e7bd971cf08055f99/xarray/coding/times.py#L32-L37

So that is why it works now.

Case closed, I guess.

I believe the issue originates in these lines:

https://github.com/pydata/xarray/blob/ac854f081d4b57d292755d3aff1476f8e2e2da11/xarray/conventions.py#L174-L175

where we multiply the num_dates with some float value and then cast to int64.

If num_dates is float32, numpy keeps float32 when multiplying with e.g. 1e9 and that somehow introduces an error. Here is a stripped version of the above:

python flat_num_dates = np.arange(100).astype('float32') n = 1e9 roundtripped = (flat_num_dates * n).astype(np.int64) / n assert np.all(flat_num_dates == roundtripped)

By the way, the factor has to be large, like 1e9. E.g. 1e6 ('ms since ...') won't give the error.

The weird thing is that the corresponding code in the current master is identical:

https://github.com/pydata/xarray/blob/50b0a69a7aa0fb7ac3afb28e7bd971cf08055f99/xarray/coding/times.py#L151-L152

I will look into why the result is still different from v0.10.0.

Also, if this really is the origin of the error, there are two easy ways to avoid this:

1. Cast flat_num_dates to float64: (flat_num_dates.astype(np.float64) * n).astype(np.int64)
2. Store _NS_PER_TIME_DELTA values as int, then numpy will do the casting.

For comparison, I added the same test to the 0.10.0 version in https://github.com/j08lue/xarray/tree/0p10-f4-time-decode, where it fails:

``` pytest -v xarray\tests\test_conventions.py

...

          self.assertArrayEqual(expected, actual_cmp)

E AssertionError: E Arrays are not equal E E (mismatch 90.0%) E x: array([datetime.datetime(2000, 1, 1, 0, 0), E datetime.datetime(2000, 1, 2, 0, 0), E datetime.datetime(2000, 1, 3, 0, 0),... E y: array(['2000-01-01T00:00:00.000000', '2000-01-02T00:00:00.003211', E '2000-01-03T00:00:00.006422', '2000-01-04T00:00:00.001245', E '2000-01-05T00:00:00.012845', '2000-01-06T00:00:00.024444',... ```

PR with a test

I'll see if I can find the time in the next few days.

@JamesSample excellent work confirming the limited scope.

I am quite sure that the issue is with this line: https://github.com/pydata/xarray/blob/v0.9.1/xarray/core/variable.py#L375-L376

dims = tuple(dim for k, dim in zip(key, self.dims) if not isinstance(k, (int, np.integer)))

When you run the test case I added in #1184, you will see that, inside __getitem__, key is a tuple with an empty array inside (key = (np.array(0, dtype=np.int64),)).

With that value, isinstance(key[0], (int, np.integer)) is False on 64 bit Windows and so dims=('time',), which has length 1.

But values = self._indexable_data[key] gives a zero-dimension array, such that values.ndim == 0...

I am unable to tell which of the two sides of the assertion expression is unexpected...

@JamesSample Thanks for restoring my credibility a bit here... But, no, I did not figure this out yet.

Appveyor apparently does not have the combination Windows 64 bit + Python 2.7 (https://ci.appveyor.com/project/jhamman/xarray-injyf/build/1.0.659), maybe that is why it does not reproduce the error?

Yes, this is Python 2. :flushed: Let's see what Travis & Co say.

@jenfly did you find a solution how to make opendap authentication work with xarray? Might be worthwhile posting it here, even though the issue has to do with the backends.

Great, safe_cast_to_index works nicely (it passes my test). I added the change to the existing PR.

Do we need to add more groupby tests to make sure the solution is safe for other cases (e.g. other data types)?

OK, here is the minimal example:

```python import xarray as xr import pandas as pd

def test_groupby_da_datetime(): """groupby with a DataArray of dtype datetime""" # create test data times = pd.date_range('2000-01-01', periods=4) foo = xr.DataArray([1,2,3,4], coords=dict(time=times), dims='time')

# create test index
dd = times.to_datetime()
reference_dates = [dd[0], dd[2]]
labels = reference_dates[0:1]*2 + reference_dates[1:2]*2
ind = xr.DataArray(labels, coords=dict(time=times), dims='time', name='reference_date')

# group foo by ind
g = foo.groupby(ind)

# check result
actual = g.sum(dim='time')
expected = xr.DataArray([3,7], coords=dict(reference_date=reference_dates), dims='reference_date')
assert actual.to_dataset(name='foo').equals(expected.to_dataset(name='foo'))

```

Making that, I found out that the problem only occurs when the DataArray used with groupby has dtype=datetime64[ns].

The problem is that we effectively feed the DataArray to pd.factorize and that goes well for most data types: Pandas checks with the function needs_i8_conversion whether it can factorize the DataArray and decides YES for our datetime64[ns]. But then in pd.factorize it fails because it tries to access DataArray.view to convert to int64.

So as I see it there are three possible solutions to this: 1. Make Pandas' pd.factorize handle our datetime DataArrays better, 2. Add an attribute .view to DataArrays, or 3. Use the solution in the above PR, which means feeding only the NumPy .values to pd.factorize.

Great to see you are pushing forward on this issue @jhamman and @shoyer. I would really have liked to contribute here, but it seems like there are quite some design choices to make, which are better left in your hands.

@jhamman sorry, I only now saw that you pointed to a previous issue on the same topic (with exactly the same considerations). I did not find that issue when I searched (for "unlimited").

You were against changing to_netcdf then. Are you still?

The above solution would not require much more than that, in set_necessary_dimensions,

def set_necessary_dimensions(self, variable): for d, l in zip(variable.dims, variable.shape): if d not in self.dimensions: self.set_dimension(d, l)

would become

def set_necessary_dimensions(self, variable): for d, l in zip(variable.dims, variable.shape): if d in self._unlimited_dimensions: l = None if d not in self.dimensions: self.set_dimension(d, l)

But maybe the encoding dict is not the way to go after all, since it contains entries per variable, while it is the dimension that must be unlimited.

Currently the dataset variables can be created in any order and their necessary dimensions created whenever needed (in the set_necessary_dimensions function). I would not like to change that logic (e.g. towards creating all dimensions required by all variables first, before adding the data variables).

So how about a new kw argument to to_netcdf, like

ds.to_netcdf(unlimited_dimensions=['time'])

or

ds.to_netcdf(dimension_unlimited={'time': True})

(the second option better for explicitly setting {'time': False})?

I think it makes sense to preserve the UNLIMITED state through read/write. In my case, I subset a netCDF file along lat and lon dimensions, leaving the time dimension untouched and would therefore expect it to pass unchanged through xarray IO (staying UNLIMITED).

However, when the dataset is indexed/subset/resampled along the unlimited dimension, it would make sense that its state is dropped. But that would require a lot of ifs and buts, so I suggest we leave that aside for now.

OK, I'd be up for taking a shot at it.

Since it is per-variable and specific to netCDF, I guess the perfect place to add this is in the encoding dictionary that you can pass to to_netcdf, right? Maybe as key unlimited? E.g.

ds.to_netcdf(encoding={'time': dict(unlimited=True)})

I need to look up whether netCDF allows for defining more than one unlimited dimension, otherwise that must throw an error.

And then it is just about passing None as length to CreateDimension, at least in netCDF4 and scipy.io.netcdf. But I did not look into how xarray handles that under the hood.

Alright, I see. In that case we should not create an entirely new syntax for this - unless there are more people than just me interested in this feature. I can just use my own little convenience function that generates a series of array indices that can be fed to isel, something like this:

python def wrap_iselkw(ds, dim, istart, istop): """Returns a kw dict for indices from `istart` to `istop` that wrap around dimension `dim`""" n = len(ds[dim]) if istart > istop: istart -= n return { dim : np.mod(np.arange(istart, istop), n) }

@jhamman You guys seem to be always ahead of me - which is great! (I just did not notice it this time before I raised the issue...)

But it would actually not be that hard to make selecting regions across the boundary of a circular dimension more convenient in xray.

It would of course be great if it worked with real coordinate (longitude) values:

ds.sel(lon=slice(-90, 15))

which, however, only works if the axis explicitly covers this range (e.g. has range -180E to 180E). If it ranges, say, from 0E to 360E, an empty DataArray is returned. There is no nice way of resolving this, I guess, because you would need to hard-code 360 or have it as a clumsy, non-CF-standard parameter and transform all values to some common range before trying to find the requested region.

But it should be less problematic with

ds.isel(lon=slice(-100, 200))

i.e. isel with negative slice start. As it is now, this also returns an empty DataArray for me (xray 0.6.0), which is in a way consistent with the sel behaviour. But asking for negative indices is actually a quite pythonic way of implying circular. So could this be implemented, now that netCDF4 supports non-contiguous indexing?

That makes sense. Great that there is an option to keep_attrs. Closing this issue.

Sorry for the confusion! The loss of attributes actually occurs when applying .mean() (rather than .load()).

See this notebook (same in nbviewer) for an example with some opendap-hosted data.

netCDF4-python uses a dimension specified by the user or an unlimited dimension it finds in the dataset. Here is the corresponding code section.

Here is a print-out of the full dataset for POP ocean model output (see that gist in nbviewer).

I can see that the heuristics exclude variables from concatenation that are associated with dimensions of other variables. But why not just exclude all that do not have a time dimension?

That is actually an excellent demonstration of the power of xray for climate data analysis, @shoyer.

Something like this (including the pandas bridge) should be included in the documentation somewhere, for example under Why xray?.

Just a thought...