I see your argument, but here's my problem. In this future where things work (assuming that NEP is accepted), and I want distributed computing with dask, units, and xarray, I have: xarray wrapping a pint array wrapping a dask array. I like composition, but that level of wrapping...feels wrong to me for some reason. Is there some elegance I'm missing here? (Other than array-like things playing together.)

And then I still need hooks in xarray so that when pint does a calculation, it can update the metadata in xarray; so it feels like we're back here anyway.

@shoyer I know elsewhere you said you weren't sure about this idea any more, but personally I'd like to push forward on this idea. Do you have problems with this approach we need to resolve? Any chance you have some preliminary code?

I think this is the right way to solve the unit issue in XArray, since at it's core unit handling is mostly a metadata operation.

We do often deal with those in my line of work as well, I just happen not to right now. But time is the one thing that already carries units, doesn't it? One can convert between various datetime64 objects and adding, subtracting, dividing timedelta64 with different units mostly works as expected (except integer division; and I haven't tried indexing with timedelta64). But I take your point about unit coordinates, and I still like the idea to provisionally add such functionality on top of an optional dependency on pint, which already has the ability to write out siunitx-latex code which then can be incorporated into plotting tools (I haven't grasped xarrays plotting tools enough yet to know how easy or difficult that part would be, though). I don't see custom dtypes with unit incorporation coming any time soon to numpy and I'm not even sure it would be the right way to go (any real dtype can have any unit; but here is not the right place to discuss that).

Good point. I didn't think of that; my coordinates happen to be either time or unitless, I think. How common is it though that the full power of a unit library is needed for coordinates? I suppose it arises with indexing, i.e. the ability to write da.sel[x=1.5 km] = value (to borrow PEP 472 syntax ;-), moreso than operations between different data arrays. With a Dataset, the coordinates would correspond to variables with their own attributes, would they not (or how else would a CF-compliant NetCDF file store units for coordinates?) so it would only require a slight expansion of the DataArray class to carry along attributes on coordinates.

When it's a bit more polished I intend to publish it somewhere, but currently several things are missing (.to(...), __rsub__, __rmul__ and friends, unit tests, some other things). I currently don't have time to add features I don't need myself (such as units on coordinates) though.

I wrote a small recipe that appears to contain basic functionality I'm looking for. There's plenty of caveats but it could be a start, if such an approach is deemed desirable at all.

``` from common import ureg # or ureg = pint.UnitRegistry()

import operator import xarray class UnitsAwareDataArray(xarray.DataArray): """Like xarray.DataArray, but transfers units """

def __array_wrap__(self, obj, context=None):
    new_var = super().__array_wrap__(obj, context)
    if self.attrs.get("units"):
        new_var.attrs["units"] = context[0](ureg(self.attrs.get("units"))).u
    return new_var

def _apply_binary_op_to_units(self, func, other, x):
    if self.attrs.get("units"):
        x.attrs["units"] = func(ureg.Quantity(1, self.attrs["units"]),
                                ureg.Quantity(1, getattr(other, "units", "1"))).u
    return x

# pow is different because resulting unit depends on argument, not on
# unit of argument (which must be unitless)
def __pow__(self, other):
    x = super().__pow__(other)
    if self.attrs.get("units"):
        x.attrs["units"] = pow(ureg.Quantity(1, self.attrs["units"]),
                               ureg.Quantity(other, getattr(other, "units", "1"))).u
    return x

for tp in ("add", "sub", "mul", "matmul", "truediv", "floordiv", "mod", "divmod"): meth = "{:s}".format(tp) def func(self, other, meth=meth, tp=tp): x = getattr(super(UnitsAwareDataArray, self), meth)(other) return self._apply_binary_op_to_units(getattr(operator, tp), other, x) func.name = meth print(func, id(func)) setattr(UnitsAwareDataArray, meth, func) del func ```

Is it not? The documentation says it's new in numpy 1.11 and we're at 1.12 now.

I tried to make a small units-aware subclass of DataArray for myself. I managed to get the right behaviour for ufuncs (I think) but somehow my subclassed _binary_op is not getting called. I guess there is some logic somewhere that leads to replacing _binary_op in a subclass doesn't work (see below). But overall, how would you feel about an optional dependency on pint with a thin layer of code in the right place?

``` class UnitsAwareDataArray(xarray.DataArray): """Like xarray.DataArray, but transfers units """

def __array_wrap__(self, obj, context=None):
    new_var = super().__array_wrap__(obj, context)
    if self.attrs.get("units"):
        new_var.attrs["units"] = context[0](ureg(self.attrs.get("units"))).u
    return new_var

@staticmethod
def _binary_op(f, reflexive=False, join=None, **ignored_kwargs):
    # NB: http://stackoverflow.com/a/26807879/974555
    x = super(UnitsAwareDataArray, UnitsAwareDataArray)(f,
        reflexive, join, **ignored_kwargs)
    # do stuff
    return x

```

Apparently __numpy_ufunc__ is too new for xarray, but it would appear that adding the right code to __array_wrap__ should work, i.e. if a units attribute is present and units are enabled through pint, evaluate something like new_var.attrs["units"] = context[0](1*ureg(self.attrs["units"])).u.

I would say using the units attribute is the most natural way to go. It could be optional and then built on top of pint, which would make it rather easy to implement:

```

ureg is a pint unit registry

y = a/b y.attrs["units"] = ureg(a.attrs["units"]) / ureg(b.attrs["units"]) ```

which if I understand the codebase correctly could be added to DataArray._binary_op. Not sure if it is similarly easy for ufuncs, is that what __numpy_ufunc__ would be for?

+1 on requiring attribute handlers to be registered at a single location because this will make for cleaner code long term.

I also prefer an approach that doesn't use context managers: I agree with @darothen's comments about the issues with their use.

Regardless of the exact implementation, I am strongly in favour of this functionality being added to xarray - I can already think of a number of very useful ways I could use it!