issue_comments: 206418244

Let me try to clarify what I mean in item 2:

Allow specification of which dims to stack.

Say you have the following dataset

``` python

ds = xr.Dataset( {'temperature': (['time','nx'], [[1,1,2,2],[2,2,3,3]] ), 'humidity': (['time','nx'], [[1,1,1,1],[1,1,1,1]] )}) ```

Now imagine you want to average humidity in temperature coordinates. (This might sound like a bizarre operation, but it is actually the foundation of a sophisticated sort of thermodynamic analysis.)

Currently this works as follows

``` python

ds = ds.set_coords('temperature') ds.humidity.groupby('temperature').sum() <xarray.DataArray 'humidity' (temperature: 3)> array([2, 4, 2]) Coordinates: * temperature (temperature) int64 1 2 3 ```

However, this sums over all time. What if you wanted to preserve the time dependence, but replace the nx coordinate with temperature. I would like to be able to say

python ds.humidity.groupby('temperature', group_over='nx').sum()

and get back a DataArray with dimensions ('time', 'temperature').

Maybe this is already possible with a sophisticated use of apply. But I don't see how to do it.