issue_comments

Thank you for the input everyone, this discussion has been very useful! Closing this issue.

@dopplershift That's a good point. It's pretty trivial to create a sympl.DataArray from an xarray.DataArray, so perhaps I should be using a decorator that will convert xarray.DataArray to sympl.DataArray whenever one is passed into a sympl call. This would be similarly easy to do in metpy. One could also write a function to convert Dataset into one that contains unit-aware DataArray objects, or an open_dataset that calls xarray.open_dataset and then does such a conversion, though I'd wonder if certain Dataset calls (e.g. mean) might undo such a conversion.

In sympl our main concerns are unit checking at the boundary of components and in properly converting units when time stepping or adding outputs of components together. Maybe sympl should only be using this DataArray subclass internally, with type conversions or wrapping when taking DataArrays into and out of its methods? That would solve a lot of our problems.

unyt may be a better choice than pint for MetPy. Like I said in Sympl we don't use pint for unit information storage, only for conversion and arithmetic, so whether it uses an ndarray subclass doesn't apply for us.

@shoyer that notation might work, thanks for pointing it out! Maybe we can think of a more natural name for the accessor ("with_units"? "keep_units"? "uarray"? "u"?). I find the "storage" of units as a string in attrs to be much cleaner than any other implementation I've seen so far (like implementations that have a unit container over an underlying array, or an array of unit-aware objects). It has the added benefit that this is how units are conventionally stored in netCDF files. I don't think it's necessary to use a class other than ndarray for data storage.

@kmpaul the main reason I stayed away from cf_units is that I had bad experiences trying to get it to build with its dependencies in the past. Particularly it's an issue that it depends on the Udunits C library, which requires MinGW to install on Windows and has generally been a headache for me. I'd much prefer a pure Python unit handling implementation. For Sympl, we don't care so much about time units, because time is stored using datetime objects (potentially from the cftime package for alternate calendars). This is also the way that time units are conventionally stored in xarray, once decoded.

It may make sense for us to use some kind of stand-alone unit-aware DataArray implementation. I'd just need to be convinced that yours is well-designed, thoroughly tested, and easy to install with pip. The main things concerning me about PhysArray are 1) As a container rather than subclass, it does not implement many of the methods of DataArray and 2) There are a few design choices I don't understand, like why calendar is always a property of a PhysArray even when it isn't storing a time, why cftime objects aren't used instead of units to manage time, and why the positive attribute is important enough for PhysArray to manage (I've never seen it in any data I've used, and it's easy to check if a DataArray is all positive or negative with a function call). We could discuss these by e-mail if you like (it's my username at uw.edu). Other possibilities are that I'll take the implementation we come up with and give it its own package, or that we'll collaborate on such a package.

For reference, here are some of the sort of methods I've been adding that aren't currently in sympl:

def multiply(self, other):
    if isinstance(other, xr.DataArray):
        result = self._dataarray * other
        result.attrs['units'] = multiply_units(self._dataarray.attrs['units'], other.attrs['units'])
    else:
        result = self._dataarray * other
        result.attrs['units'] = self._dataarray.attrs['units']
    return result

def divide(self, other):
    print(self._dataarray, other)
    if isinstance(other, xr.DataArray):
        result = self._dataarray / other
        result.attrs['units'] = divide_units(self._dataarray.attrs['units'], other.attrs['units'])
    else:
        result = self._dataarray / other
        result.attrs['units'] = self._dataarray.attrs['units']
    return result

def add(self, other):
    result = self._dataarray + other.sympl.to_units(self._dataarray.attrs['units'])
    result.attrs['units'] = self._dataarray.attrs['units']
    return result

def subtract(self, other):
    result = self._dataarray - other.sympl.to_units(self._dataarray.attrs['units'])
    result.attrs['units'] = self._dataarray.attrs['units']
    return result