issue_comments

My personal hope is to keep this thread open just for the record. But given the non-activity on the numpy end, I honestly can't promise any resolution to this issue in the near future. Thanks!

PS I persist because some people do seem to appreciate that PR and have forked it for their own use :)

So, as per suggestion by Stephan, I submitted the PR to numpy instead: https://github.com/numpy/numpy/pull/9211

Despite having possibly the highest number of comments of all active numpy PRs and passing all the tests, it's been sitting in limbo for the last few months. There seems to be quite some uncertainty about what to do with it, and the discussion has gone off tangent a bit to more upstream issues.

My hope is still to have it eventually merged into numpy, so that it could be easily ported to xarray. It's proven to be quite useful to my coworkers, as well as some of the numpy users. I believe it'll also serve xarray well. Thank you all for your patience.

This is getting close to the fundamental paradigm of xarray. Would it entail some massive code migration (and user protests) if we change how we express da.coords['x'](by removing point coordinates)?

I'm thinking of a new API that compares only the values and indexing coordinates between two objects, but struggle to come up with a simple, descriptive name........

I'm just having a hard time wrapping my head around the definitions here. da['some_attr'] = 0 sets a point coordinate, and this point coordinate shows up in da.coords, but not da.dims, because da.dims only pertains to indexing coordinates. Why is this point coordinate attached to every other coordinate then? In other words, why would da.coords['x'] contain some_attr as well?

Somehow it just doesn't feel right intuitively for some point dimension to mess up my comparison between two 1-D coordinates, but I might be missing some important use cases here.

Thank you very much for your patience!

BTW, I think I messed up my co-worker's name in whats-new.rst. The spelling is right (Kyle Heuton), but the extra indentation shouldn't be there. As a result his name appears a bit out of place. Would be very much obliged if you could fix that in your next push to master. Gracias.

api.rst should've been a checklist item :)

If that's the case then I have nothing more to add. I've run flake8 on the individual files I touched I as well, such as flake8 xarray/core/dataset.py, and none of my changes are flagged.

git diff upstream/master | flake8 --diff only complains about a bunch of lines in whats-new.rst and reshaping.rst. But these complaints don't make any sense to me. Some examples: doc/reshaping.rst:195:1: E112 expected an indented block doc/reshaping.rst:203:1: E101 indentation contains mixed spaces and tabs doc/reshaping.rst:203:1: W191 indentation contains tabs doc/reshaping.rst:203:2: E113 unexpected indentation I've played around with it a bit but not seen any changes to the complaints. What is going on? Should I even worry about the .rst files?

Haha, brilliant.

I just made some futile attempts to find that checklist, to no avail. Mind helping me with a link? :)

Today is the first time I've heard of flake8. Is this how you guys standardize quality checks? If so, do you run once for python 3.5 and once for 2.7?

Looks ready now.

Think I just copied something over from variable.py and changed the message a little bit...

From numpy 1.11.0 release notes: ``` np.lexsort now supports arrays with object data-type

The function now internally calls the generic npy_amergesort when the type does not implement a merge-sort kind of argsort method. ```

1.10.4 failed with the merge sort error. 1.11.0 failed with some error that had nothing to do with merge sort (Reason: image not found. Make sure your test modules/packages have valid Python names.) Starting from 1.11.1 it works fine.

What's the best course of action now?

If I downgrade to numpy 1.10, I get that error. Once upgraded to 1.12, it goes away......

Perhaps it's finally ready for prime time?

I've added a tiny bit of extra docstring just to ease my discomfort about the NaN sort. If you feel it's redundant I can remove it.

Almost time to merge this baby in? :)

This has impacted our workflow as well. Please let us know if we can help.

Weird. Travis seems to fail in places unrelated to my changes.

I'm rather agnostic to sortby vs. sort_by.

Sounds good. As I'm writing the type-checking code I run into this question: Why would I have a xarray.core.dataset.DataVariables object as input? A DataVariables object could contain multiple DataArrays, which makes the interpretation a bit unclear. In my mind it should only be 1.) name(s) of existing index coords, or 2.) 1D DataArray(s) with dim in self.dims

Just did. Got it to work. Muchas Gracias.

What would the signature of sortby() be then? On our end we just want a more intuitive way to sort by dimension labels, so now I have sort_index(self, dims, ascending=True). sortby(), based on your description, seems like a separate method. Or any suggestion on how we can marry the two into something coherent?

A couple of things:

1.) Upon a little thinking I believe sort_values() doesn't make much sense, so I'm only working on sort_index()'. 2.) the way I handle theinplacekwarg is by ``` if inplace: self = self.isel(**{d: self.indexes[d].argsort() if ascending else self.indexes[d].argsort()[::-1] for d in dimensions}) else: return self.isel(**{d: self.indexes[d].argsort() if ascending else self.indexes[d].argsort()[::-1] for d in dimensions}) ``` But when I run ``` ds.sort_index(dims=['x', 'y'], inplace=True) ``` nothing changes. If I put apdb.set_trace()` right below the self = self*** I can evaluate self and see that it's what I want it to be. But somehow that assignment is not realized to the higher level. Any quick pointer?

When I do conda install netcdf4, it ends up downgrading my numpy from 1.12.0 to 1.11.3, which triggers other dependency-related problems in my code. Is there a different solution to this problem?

On our end, we currently do the following when we need to sort by axis label (lat/lon in this case): da.reindex(indexers={'lat':sorted(da.coords['lat'].values), 'lon':sorted(da.coords['lon'].values)}) Upon first glance of the source code I think our approach goes down different code path than your .isel() approach. The most obvious difference, from a user's stand point, is probably that .reindex() returns a new object, whereas .isel() returns a view (typically). In Pandas, both sort_index() and sort_values() seem to return new objects.

We'd be happy to contribute to an xarray version of sort_index() and sort_values(). The first question is, which one would be the more robust and computationally efficient code path to take?

Cool, let me toss it out on the numpy board and see what they think.

This functionality is quite important to us as well. How might we help? Should I simply fork groupby-multiple-args?

tears in my eyes.......

checking in on this PR...

time to merge this guy?

I'm not sure why continuous integration check is taking forever.....

So I addressed some of the comments, but not all. Please let me know what you think.

I'm curious, if you put a from .computation import apply_ufunc inside ops.py, would you not get some circular ImportError? Seems difficult to get out of this circle......

Ok I'll give it a shot. Will touch base when I run into roadblocks.

So I took at quick look at the commits in #964. It's not entirely clear to me how one can easily add a join argument to fillna. Should I keep my current commits, and just submit a PR to master once #964 is merged in, and then we could see how it goes from there?

Checking in on how to move forward from here...I feel it's pretty close to a PR...

Currently, ds0.combine_first(ds1) gives exactly the same result as xr.merge([ds0, ds1]). But for data arrays it still offers something new.

Either 1.) my combine_first should be doing something different with datasets, or 2.) we don't need a combine_first for datasets, or 3.) change xr.merge so that when applied to data arrays, it creates a new data array with outer-join and fillna (by chaining combine_first down the list of data arrays).

I suppose I can save one line by getting rid of the duplicate f = _func_slash_method_wrapper(fillna), but besides that I'm not sure what's the best way to refactor this.

If the behaviors check out, then this branch might be ready for a PR.

So these are my _fillna and _combine_first in ops.inject_binary_ops:

``` f = _func_slash_method_wrapper(fillna) method = cls._binary_op(f, join='left', fillna=True) setattr(cls, '_fillna', method)

f = _func_slash_method_wrapper(fillna) method = cls._binary_op(f, join='outer', fillna=True) setattr(cls, '_combine_first', method) ```

Within dataarray.py and dataset.py, combine_first(self, other) simply returns self._combine_first(other). This code path renders the test results you saw earlier.

Given this construct, I'm not sure how to do the refactor like you mentioned. Perhaps a few more pointers to the right direction? :)

hmm, so what would be an expected behavior of ds.combine_first?

If I have

```

ds0 <xarray.Dataset> Dimensions: (x: 2, y: 2) Coordinates: * x (x) |S1 'a' 'b' * y (y) int64 -1 0 Data variables: ds0 (x, y) int64 0 0 0 0 ds1 <xarray.Dataset> Dimensions: (x: 2, y: 2) Coordinates: * x (x) |S1 'b' 'c' * y (y) int64 0 1 Data variables: ds1 (x, y) int64 1 1 1 1 ```

I get

```

ds0.combine_first(ds1) <xarray.Dataset> Dimensions: (x: 3, y: 3) Coordinates: * x (x) object 'a' 'b' 'c' * y (y) int64 -1 0 1 Data variables: ds0 (x, y) float64 0.0 0.0 nan 0.0 0.0 nan nan nan nan ds1 (x, y) float64 nan nan nan nan 1.0 1.0 nan 1.0 1.0 ```

and changing the order to ds1.combine_first(ds0) just flips the order of the data_vars, but the cell values of the data_vars remain the same.

This is done essentially by adding a _combine_first to ops.py that mimics _fillna, except join='outer'.

woohoo!

Any suggestion on improvement?

Alright, if that's desired, I'd just change code in apply_over_both() of dataset._calculate_binary_op to react to different kwarg values, and adding a method in test_dataset.py. Should be able to push before the weekend is over.

Just checking in. Is there something else that needs to be added to this PR before merge?

So I spent some time looking at dataset._calculate_binary_op but couldn't quite come up with what I wanted. After banging my head against the wall a bit this is what I have:

def combine_first(left, right): la, ra = xr.align(left, right, join='outer', copy=False) # should copy=True? la, ra = la.where(ra.isnull() | ra.notnull()), ra.where(la.isnull() | la.notnull()) ra.values[la.notnull().values] = la.values[la.notnull().values] return ra

And it seems to work. My test cases are

```

l_2d <xarray.DataArray (x: 2, y: 2)> array([[1, 1], [1, 1]]) Coordinates: * x (x) |S1 'a' 'b' * y (y) int64 -2 0 r_2d <xarray.DataArray (x: 2, y: 2)> array([[0, 0], [0, 0]]) Coordinates: * x (x) |S1 'b' 'c' * y (y) int64 0 2 ar_1d <xarray.DataArray (x: 3)> array([4, 5, 6]) Coordinates: * x (x) |S1 'a' 'b' 'd' ```

and here are the results:

```

combine_first(l_2d, r_2d) <xarray.DataArray (x: 3, y: 3)> array([[ 1., 1., nan], [ 1., 1., 0.], [ nan, 0., 0.]]) Coordinates: * x (x) object 'a' 'b' 'c' * y (y) int64 -2 0 2 combine_first(r_2d, l_2d) <xarray.DataArray (x: 3, y: 3)> array([[ 1., 1., nan], [ 1., 0., 0.], [ nan, 0., 0.]]) Coordinates: * x (x) object 'a' 'b' 'c' * y (y) int64 -2 0 2 combine_first(l_2d, ar_1d) <xarray.DataArray (x: 3, y: 2)> array([[ 1., 1.], [ 1., 1.], [ 6., 6.]]) Coordinates: * x (x) object 'a' 'b' 'd' * y (y) int64 -2 0 combine_first(ar_1d, l_2d) <xarray.DataArray (x: 3, y: 2)> array([[ 4., 4.], [ 5., 5.], [ 6., 6.]]) Coordinates: * x (x) object 'a' 'b' 'd' * y (y) int64 -2 0 ```

I don't like the fact that I have to access .values, and the use of .where is slightly wonky. But this is definitely the cleanest working solution I have thus far. Any suggestion to improve is welcome.

Thanks for the reply. Would it make things easier if c, d = xr.align(a, b, join="outer") explicitly broadcasts missing dimensions in all returned values? Currently, if b is missing a dimension, then d would also miss a dimension:

```

a <xarray.DataArray (x: 2, y: 2)> array([[1, 1], [1, 1]]) Coordinates: * x (x) |S1 'a' 'b' * y (y) int64 -2 0 b <xarray.DataArray (x: 3)> array([4, 5, 6]) Coordinates: * x (x) |S1 'a' 'b' 'd' c, d = xr.align(a, b, join="outer") c <xarray.DataArray (x: 3, y: 2)> array([[ 1., 1.], [ 1., 1.], [ nan, nan]]) Coordinates: * x (x) object 'a' 'b' 'd' * y (y) int64 -2 0 d <xarray.DataArray (x: 3)> array([4, 5, 6]) Coordinates: * x (x) object 'a' 'b' 'd' ```

Another test:

```

left <xarray.DataArray (x: 2, y: 2)> array([[1, 1], [1, 1]]) Coordinates: * x (x) |S1 'a' 'b' * y (y) int64 -2 0 right <xarray.DataArray (x: 2, y: 2)> array([[0, 0], [0, 0]]) Coordinates: * x (x) |S1 'b' 'c' * y (y) int64 0 2 combine_first(left, right) <xarray.DataArray (x: 3, y: 3)> array([[ 1., 1., nan], [ 1., 1., 0.], [ nan, 0., 0.]]) Coordinates: * x (x) object 'a' 'b' 'c' * y (y) int64 -2 0 2 combine_first(right, left) <xarray.DataArray (x: 3, y: 3)> array([[ 1., 1., nan], [ 1., 0., 0.], [ nan, 0., 0.]]) Coordinates: * x (x) object 'a' 'b' 'c' * y (y) int64 -2 0 2 ```

Here's my somewhat hacky stab at it, not optimized for speed/efficiency:

def combine_first(left, right): """ Takes 2 data arrays, performs an outer-concat, using values from the right array to fill in for missing coordinates in the left array. """ l_aligned, r_aligned = xr.align(left, right, join="outer") temp = l_aligned + r_aligned # hack temp.values[temp.notnull().values] = np.nan # now template is all nan # insert non-nan values from right array first temp.values[r_aligned.notnull().values] = r_aligned.values[r_aligned.notnull().values] # insert values from left array, overwriting those from right array temp.values[l_aligned.notnull().values] = l_aligned.values[l_aligned.notnull().values] return temp

And the result:

```

ar1 = xr.DataArray([4,5,6],[('x',[1,2,3])]) ar2 = xr.DataArray([[7,8,9],[10,11,12],[13,14,15]],[('x',[1,12,13]),('y',[0,5,6])]) ar1 <xarray.DataArray (x: 3)> array([4, 5, 6]) Coordinates: * x (x) int64 1 2 3 ar2 <xarray.DataArray (x: 3, y: 3)> array([[ 7, 8, 9], [10, 11, 12], [13, 14, 15]]) Coordinates: * x (x) int64 1 12 13 * y (y) int64 0 5 6 temp = combine_first(ar1, ar2) temp <xarray.DataArray (x: 5, y: 3)> array([[ 4., 5., 6.], [ 4., 5., 6.], [ 4., 5., 6.], [ 10., 11., 12.], [ 13., 14., 15.]]) Coordinates: * x (x) int64 1 2 3 12 13 * y (y) int64 0 5 6 temp = combine_first(ar2, ar1) temp <xarray.DataArray (x: 5, y: 3)> array([[ 7., 8., 9.], [ 4., 5., 6.], [ 4., 5., 6.], [ 10., 11., 12.], [ 13., 14., 15.]]) Coordinates: * x (x) int64 1 2 3 12 13 * y (y) int64 0 5 6 ```

Would this be the behavior we want from xarray's combine_first? Not entirely analogous to the Pandas function. Maybe rename it?

So I fooled around with Pandas' combine_first:

```

df1 = pd.DataFrame({'x':[1,2,3],'z':[4,5,6]}).set_index('x') df1
z x
1 4 2 5 3 6 df2 = pd.DataFrame({'x':[1,12,13],'y':[0,5,6],'z':[7,8,9]}).set_index(['x','y']) df2 z x y
1 0 7 12 5 8 13 6 9 df1.combine_first(df2) z x y
1 0 4.0 12 5 8.0 13 6 9.0 ```

and was surprised that the indexes were not "outer-joined". Is this the behavior xarray wants to emulate?

As a mockup for xr.combine_first(arr1, arr2), I was thinking about using align([arr1, arr2], join="outer") to set up the template, and then go into the template to set the values right. Is that a sound approach? Muchas gracias.

Is combine_first already implemented? I can't find it in the source code, neither could I find ops.choose.

Out of curiosity, how do I test inner-join vs. outer-join using only 1-D arrays? I can certainly simplify the current tests to be 2D + 1D.

Yup, did that. Thank you. Will request to have you as a reviewer. Much obliged.

Thanks. I will soon submit a pull request from my fork. Which branch should it go to? It seems that many people submit pull requests to master here.

Let me exhaust a few other ideas first. I'll definitely share my thoughts here first before making any commit. Thanks.

My opinion is that the nan has got to go. If we want to (1) maintain pandas-consistency and (2) use bottleneck without mucking it up, then I think we need to add some logic in either rolling.reduce() or rolling._center_result().

So here's my failed attempt:

``` def reverse_and_roll_1d(data, window_size, min_periods=1): """ Implements a concept to fix the end-of-array problem with xarray.core.rolling._center_shift(), by 1.) take slice of the back-end of the array 2.) flip it 3.) compute centered-window arithmetic 4.) flip it again 5.) replace back-end of default result with (4)

:param DataArray data: 1-D data array, with dim name 'x'.
:param int window_size: size of window.
"""
# first the default way to computing centered window
r = data.rolling(x=window_size, center=True, min_periods=min_periods)
avg = r.mean()
# now we need to fix the back-end of the array
rev_start = len(data.x) # an index
rev_end = len(data.x) - window_size - 1 \
                 if len(data.data) > window_size \
                else None  # another index
tail_slice = slice(rev_start, rev_end, -1) # back end of array, flipped
r2 = data[dict(x=tail_slice)].\
    rolling(x=window_size, center=True, min_periods=min_periods)
avg[dict(x=slice(-window_size+1, None))] = \
    r2.mean()[dict(x=slice(window_size-2, None, -1))] # replacement

return avg

```

This algorithm is consistently 8 times slower than pd.DataFrame.rolling(), for various 1d array sizes.

I'm open to ideas as well :)