Repository: xgcm/xhistogram

Issue:
axes not used for making hist got shifted
Say we have two xarray.DataArrays with dims (t:4, z:10, x:100, y:100). If we make histagram based on dimensions x and y, what we have in the result is  (t:10, z:4, x_bin, y_bin).  The t and z get exchanged. 


```python
import numpy as np
import xarray as xr
from xhistogram.xarray import histogram


# create data for testing

time_axis = range(4)
depth_axis = range(10)
X_axis = range(30)
Y_axis = range(30)

dat1 = np.random.randint(low = 0, high = 100,size=(len(time_axis),len(depth_axis),len(X_axis),len(Y_axis)))
array1 = xr.DataArray(dat1, coords = [time_axis,depth_axis,X_axis,Y_axis], dims = ['time','depth','X','Y'])

dat2 = np.random.randint(low = 0, high = 50,size=(len(time_axis),len(depth_axis),len(X_axis),len(Y_axis)))
array2 = xr.DataArray(dat2, coords = [time_axis,depth_axis,X_axis,Y_axis], dims = ['time','depth','X','Y'])

# create bins and rename arrays

bins1 = np.linspace(0, 100, 50)
bins2 = np.linspace(0,50,25)

array1 = array1.rename('one')
array2 = array2.rename('two')

result= histogram(array1,array2,dim = ['X','Y'] , bins=[bins1,bins2])
```

The dimensions of result is (time: 10, depth: 4, one_bin: 49, two_bin: 24) instead of (time: 4, depth: 10, one_bin: 49, two_bin: 24). Is this a bug of the code or just my misusing of the function?

Discussion:

Patch:
diff --git a/xhistogram/xarray.py b/xhistogram/xarray.py
index 230a44c..c4d41d8 100644
--- a/xhistogram/xarray.py
+++ b/xhistogram/xarray.py
@@ -4,6 +4,7 @@ Xarray API for xhistogram.
 
 import xarray as xr
 import numpy as np
+from collections import OrderedDict
 from .core import histogram as _histogram
 
 
@@ -95,11 +96,11 @@ def histogram(*args, bins=None, dim=None, weights=None, density=False,
 
     # roll our own broadcasting
     # now manually expand the arrays
-    all_dims = set([d for a in args for d in a.dims])
-    all_dims_ordered = list(all_dims)
+    all_dims = [d for a in args for d in a.dims]
+    all_dims_ordered = list(OrderedDict.fromkeys(all_dims))
     args_expanded = []
     for a in args:
-        expand_keys = all_dims - set(a.dims)
+        expand_keys = [d for d in all_dims_ordered if d not in a.dims]
         a_expanded = a.expand_dims({k: 1 for k in expand_keys})
         args_expanded.append(a_expanded)
 
@@ -118,7 +119,7 @@ def histogram(*args, bins=None, dim=None, weights=None, density=False,
         weights_data = None
 
     if dim is not None:
-        dims_to_keep = [d for d in a_dims if d not in dim]
+        dims_to_keep = [d for d in all_dims_ordered if d not in dim]
         axis = [args_transposed[0].get_axis_num(d) for d in dim]
     else:
         dims_to_keep = []
@@ -129,11 +130,19 @@ def histogram(*args, bins=None, dim=None, weights=None, density=False,
 
     # create output dims
     new_dims = [a.name + bin_dim_suffix for a in args[:N_args]]
-    bin_centers = [0.5*(bin[:-1] + bin[1:]) for bin in bins]
-    coords = {name: ((name,), bin_center, a.attrs)
-              for name, bin_center, a in zip(new_dims, bin_centers, args)}
     output_dims = dims_to_keep + new_dims
 
+    # create new coords
+    bin_centers = [0.5*(bin[:-1] + bin[1:]) for bin in bins]
+    new_coords = {name: ((name,), bin_center, a.attrs)
+                  for name, bin_center, a in zip(new_dims, bin_centers, args)}
+
+    old_coords = {name: a0[name]
+                  for name in dims_to_keep if name in a0.coords}
+    all_coords = {}
+    all_coords.update(old_coords)
+    all_coords.update(new_coords)
+
     # CF conventions tell us how to specify cell boundaries
     # http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#cell-boundaries
     # However, they require introduction of an additional dimension.
@@ -143,7 +152,8 @@ def histogram(*args, bins=None, dim=None, weights=None, density=False,
                   for name, bin_edge, a in zip(edge_dims, bins, args)}
 
     output_name = '_'.join(['histogram'] + [a.name for a in args[:N_args]])
-    da_out = xr.DataArray(h_data, dims=output_dims, coords=coords,
+
+    da_out = xr.DataArray(h_data, dims=output_dims, coords=all_coords,
                           name=output_name)
     return da_out
 


Test patch:
diff --git a/xhistogram/test/test_xarray.py b/xhistogram/test/test_xarray.py
index 9b46dd4..0018f0d 100644
--- a/xhistogram/test/test_xarray.py
+++ b/xhistogram/test/test_xarray.py
@@ -94,3 +94,31 @@ def test_weights(ones, ndims):
                 for d in combinations(dims, nc+1):
                     h = histogram(ones, weights=weights, bins=[bins], dim=d)
                     _check_result(h, d)
+
+
+# test for issue #5
+def test_dims_and_coords():
+    time_axis = np.arange(4)
+    depth_axis = np.arange(10)
+    X_axis = np.arange(30)
+    Y_axis = np.arange(30)
+
+    dat1 = np.random.randint(low=0, high=100,
+                             size=(len(time_axis), len(depth_axis),
+                                   len(X_axis), len(Y_axis)))
+    array1 = xr.DataArray(dat1, coords=[time_axis,depth_axis,X_axis,Y_axis],
+                          dims=['time', 'depth', 'X', 'Y'], name='one')
+
+    dat2 = np.random.randint(low=0, high=50,
+                             size=(len(time_axis), len(depth_axis),
+                                   len(X_axis), len(Y_axis)))
+    array2 = xr.DataArray(dat2, coords=[time_axis,depth_axis,X_axis,Y_axis],
+                          dims=['time','depth','X','Y'], name='two')
+
+    bins1 = np.linspace(0, 100, 50)
+    bins2 = np.linspace(0,50,25)
+
+    result = histogram(array1,array2,dim = ['X','Y'] , bins=[bins1,bins2])
+    assert result.dims == ('time', 'depth', 'one_bin', 'two_bin')
+    assert result.time.identical(array1.time)
+    assert result.depth.identical(array2.depth)