awesome work @keewis !

Phew thanks a lot @keewis!

Great! Quite an effort @keewis —thank you!

I would think so, too. Go ahead.

Edit: seems like other PRs have the same problem, see #3409

I think these tests are unrelated and we can merge?

```python

=================================== FAILURES =================================== ____ testdatetime_reduce[False] ______

dask = False

@arm_xfail
@pytest.mark.parametrize("dask", [False, True])
def test_datetime_reduce(dask):
    time = np.array(pd.date_range("15/12/1999", periods=11))
    time[8:11] = np.nan
    da = DataArray(np.linspace(0, 365, num=11), dims="time", coords={"time": time})

    if dask and has_dask:
        chunks = {"time": 5}
        da = da.chunk(chunks)

    actual = da["time"].mean()

  assert not pd.isnull(actual)

E AssertionError: assert not True E + where True = <function isna at 0x7fcbdefdbae8>(<xarray.DataArray 'time' ()>\narray('NaT', dtype='datetime64[ns]')) E + where <function isna at 0x7fcbdefdbae8> = pd.isnull

xarray/tests/test_duck_array_ops.py:288: AssertionError ____ test_datetime_reduce[True] ______

dask = True

@arm_xfail
@pytest.mark.parametrize("dask", [False, True])
def test_datetime_reduce(dask):
    time = np.array(pd.date_range("15/12/1999", periods=11))
    time[8:11] = np.nan
    da = DataArray(np.linspace(0, 365, num=11), dims="time", coords={"time": time})

    if dask and has_dask:
        chunks = {"time": 5}
        da = da.chunk(chunks)

    actual = da["time"].mean()

  assert not pd.isnull(actual)

E AssertionError: assert not True E + where True = <function isna at 0x7fcbdefdbae8>(<xarray.DataArray 'time' ()>\narray('NaT', dtype='datetime64[ns]')) E + where <function isna at 0x7fcbdefdbae8> = pd.isnull ```

done

Also please move the what's new blob from 0.14.0 to 0.14.1

whatsnew.rst: ``` - Added integration tests against pint https://pint.readthedocs.io/_.

.. note::

At the moment of writing, these tests *as well as the ability to use pint in general*
require `a highly experimental version of pint
<https://github.com/andrewgsavage/pint/pull/6>`_ (install with 
``pip install git+https://github.com/andrewgsavage/pint.git@refs/pull/6/head)``.
Even with it, interaction with non-numpy array libraries, e.g. dask or sparse, is broken.

installing.rst - pint <https://pint.readthedocs.io/>_: for units of measure.

.. note::

At the moment of writing, xarray requires a `highly experimental version of pint 
<https://github.com/andrewgsavage/pint/pull/6>`_ (install with
``pip install git+https://github.com/andrewgsavage/pint.git@refs/pull/6/head)``.
Even with it, interaction with non-numpy array libraries, e.g. dask or sparse, is broken.

```

I don't really know what to write about

readers will be much more interested in what they can actually do with pint

because a feature list seems to be a bit too long to me.

re installing.rst: is adding (experimental) enough, or should the line with pip also appear there?

@keewis looks good. Also something equivalent in installing.rst on the line that mentions pint, please

I have been using a PR to hgrecco/pint#764 because that has been the version that had the most tests pass. Is that okay, or should I use the original PR instead?

So something like this?

Added integration tests against pint <https://pint.readthedocs.io/>_. Note that these tests currently require a special version of pint <https://github.com/andrewgsavage/pint/pull/6> (install with pip install git+https://github.com/andrewgsavage/pint.git@refs/pull/6/head) to run and are otherwise skipped. Also note that at the moment of writing, there are issues when pint interacts with non-numpy array libraries, e.g. dask or sparse.

Edit: this one works

Added integration tests against pint <https://pint.readthedocs.io/>_. Note that these tests currently require a special version of pint to run and are otherwise skipped. Also note that at the moment of writing, there are issues when pint interacts with non-numpy array libraries, e.g. dask or sparse.

(remove emphasis in actual version)

readers will be much more interested in what they can actually do with pint

👍 , let's add a sentence

@keewis maybe you could clarify where to find this special version and how to get it (pip install https+git://...)? Also I think readers will be much more interested in what they can actually do with pint and what they can't, as opposed to what CI tests we're running.

what do you think of this for whats-new.rst?

Added integration tests against pint <https://pint.readthedocs.io/>_. Note that these tests currently require a special version of pint to run. Also note that at the moment of writing, there are issues when pint interacts with non-numpy array libraries, e.g. dask or sparse.

@keewis could you tweak whatsnew.rst and installing.rst to set the expectations (very low) for people using the current latest available release of pint?

OK great! @crusaderky any final thoughts? Anyone else?

I assumed this to be the same as datetime, which means object labels instead of strings. So now only the top-level functions and dataset are missing.

👍

I updated the PR so that all the DataArray methods I found in the docs have a test, the only exception is to_unstacked_dataset where I don't know what the result should look like: Consider ```python

data = [[0, 10, 20], [30, 40, 50]] x = [0, 1] y = [0, 1, 2] arr = xr.DataArray(data=data, coords={"x": x, "y": y}, dims=("x", "y")) arr.stack(z=("x", "y")).to_unstacked_dataset(dim="z") <xarray.Dataset> Dimensions: (y: 3) Coordinates: * y (y) int64 0 1 2 Data variables: 0 (y) int64 0 10 20 1 (y) int64 30 40 50 `` and supposing that indexes keep the units, what should the result of the above, but with units onx`, look like?

Edit: I assumed this to be the same as datetime, which means object labels instead of strings. So now only the top-level functions and dataset are missing.

I'd like to leave these two and replacing assert_equal_with_units with either xarray.testing.assert_identical or xarray.testing.assert_allclose or a combination of both to new PRs, so now this should be ready for review.

Fixed CI, What's New, and Installing doc

that's fine with me. I resolved the merge conflicts (along with other stuff), so only whats-new.rst and sorting out what is to be done in new PRs are left to do.

@keewis - I'm wondering if it would be worth merging this in its current state and leaving the rest of the checklist for a follow up PR? No harm in working incrementally on this. We just need to add a note in whats-new.rst and resolve a few merge conflicts. Thoughts?

In general we've tried to make identical as exactly identical, including types etc, so I think different units would return False there

I'm currently working on the tests for the comparison methods (DataArray.equals, DataArray.identical and DataArray.broadcast_equals) and I'm trying to figure out what equals and identical return. equals is easy: just make sure the wrapped quantities are equal. identical is a stricter version of equals, so I was wondering whether it should also return False if the values are equal, but in compatible units (like 10 cm and 0.1 m)? Any thoughts on this?

@keewis Thank you for catching this! I've pushed an update to https://github.com/andrewgsavage/pint/pull/6 with your suggested changes. Raising an error on incompatible/missing units is definitely something that will be need to be added, but it may take some re-thinking the current implementations and how convert_to_consistent_units is used. My instinct is to auto-wrap non-Quantities (of known types on the type casting hierarchy) as dimensionless (as in https://github.com/hgrecco/pint/pull/764#issuecomment-523272038), and let pint's usual unit comparison checks handle the rest. But, knowing exactly what to do here and when may be something that would have to wait until some other discussions with pint are cleared up a little more (like https://github.com/hgrecco/pint/pull/764 vs. https://github.com/hgrecco/pint/pull/875 and https://github.com/hgrecco/pint/issues/845 / https://github.com/hgrecco/pint/issues/878)?

I'm thinking that discussion regarding bugs with pint's upcoming __array_function__ may fit best on https://github.com/andrewgsavage/pint/pull/6 right now, since https://github.com/hgrecco/pint/pull/764 seems to have halted in favor of https://github.com/hgrecco/pint/pull/875? (That is, until a new PR is made that follows up on https://github.com/hgrecco/pint/pull/875 to replicate what was done in https://github.com/hgrecco/pint/pull/764 and https://github.com/andrewgsavage/pint/pull/6.)

I'm not sure where to put this, so I'll post this here.

While working on the numpy methods part and using any of the pint versions implement __array_function__ (hgrecco/pint#764 and andrewgsavage/pint#6) , I encountered the following: ```python

array = np.array([1.3, 3.1, 7.5]) array.searchsorted(2.4) 1 (array * ureg.m).searchsorted(2.4 * ureg.m) Traceback (most recent call last): File "<stdin>", line 1, in <module> File ".../pint/quantity.py", line 1548, in searchsorted raise DimensionalityError('dimensionless', self._units) pint.errors.DimensionalityError: Cannot convert from 'dimensionless' to 'meter' (array * ureg.m).searchsorted(2.4) <Traceback> pint.errors.DimensionalityError: Cannot convert from 'dimensionless' to 'meter' (array * ureg.m).searchsorted(2.4 * ureg.s) <Traceback> pint.errors.DimensionalityError: Cannot convert from 'dimensionless' to 'meter' I believe that is because in `BaseQuantity.searchsorted`, `self.__class__` is used to check the type, which may be either `QuantitySequence` or `QuantityScalar`. If `self` is a sequence and a scalar is passed to the method, this fails. To fix this, replacing `self.__class__` with `BaseQuantity` seems to work.python np.searchsorted(array * ureg.m, 2.4 * ureg.m) <Quantity(1, 'meter')> np.searchsorted(array * ureg.m, 2.4) <Quantity(1, 'meter')> np.searchsorted(array * ureg.m, 2.4 * ureg.s) <Traceback> pint.errors.DimensionalityError: Cannot convert from 'second' ([time]) to 'meter' ([length]) `` I don't know what to do here (using the code from andrewgsavage/pint#6 I tried movingsearchsortedinto the same category asisclose` -- this seems to work partially, but passing non-quantities is still allowed), but the result should be indices and incompatible units (or no unit at all) should raise an error.

Both the method and the function work fine with current master.

me hitting that issue is probably a bug in how xr.Variable handles np.full_like (the unit gets stripped), so no worries.

@keewis Thank you for pointing that out, I forgot to mention that right now mixed types are not handled by https://github.com/hgrecco/pint/pull/764 (see https://github.com/hgrecco/pint/pull/764#issuecomment-523272038).

I'm not sure whether that belongs into hgrecco/pint#882, so I'll ask this here: In the case of full_like with other being a ndarray and fill_value having a unit, numpy complains about np.copyto not being implemented for this type of operation: ```python

np.full_like(np.arange(10), fill_value=1 * ureg.m) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "<array_function internals>", line 6, in full_like File ".../lib/python3.7/site-packages/numpy/core/numeric.py", line 401, in full_like multiarray.copyto(res, fill_value, casting='unsafe') File "<array_function internals>", line 6, in copyto TypeError: no implementation found for 'numpy.copyto' on types that implement array_function: [<class 'numpy.ndarray'>, <class 'pint.quantity.build_quantity_class.\<locals>.QuantityScalar'>]

yeah, well, I can create that index outside of pytest with python pd.Index(list(np.arange(10) * ureg.m)) so the performance is that of an index on a python list. How that happens is beyond me, though.

@keewis My inclination is to think of the units as part of the data, and that, for example, zeros_like returning an array of bare zeros instead of zeros-with-units is reasonably intuitive. But since this is likely not everyone's perspective, I think raising an issue on pint's end would be good.

I'm not sure about the indexing behavior. From the sounds of your prior comment, it works because it functions as an "object-type" index? If so, it may cause a decent hit on performance. I'd definitely want to hear others' thoughts on it too.

About the indexing: this might have been accidental (or a bug in pandas) and is probably a bad idea, but do you think it might be worth to somehow actively trigger this behavior in pd.Index for pint arrays (and possibly others) until the index refactoring from https://github.com/pydata/xarray/issues/1603#issuecomment-511126208 is done?

I like your approach because it avoids all those issues, but stripping the units might be a surprise. Then again, what does a function named *_like promise? From the docs, it seems that only shape and dtype are guaranteed, so the question is whether or not the unit is part of the data or part of the type?

@keewis In https://github.com/andrewgsavage/pint/pull/6, I implemented zeros_like, ones_like and empty_like to return base ndarrays rather than quantities, and full_like to have the unit of fill_value. This seemed like the most sensible behavior to me, since it avoids the ambiguities you mention in the full_like case if they were based on the unit of other, and for many of my own use cases, I've often wanted a zeros/ones array with a different unit than the quantity whose shape I am basing it off of.

Does this behavior seem reasonable to you? Also, would this be something that should be cleared up with an issue on pint's end?

I tried adding tests for zeros_like, ones_like and full_like, but I'm not sure what the expected result should look like, at least for the last one. I would assume zeros_like and ones_like return a quantity with the same unit as other, but what is the unit for full_like? The unit of other, or the unit of fill_value? Also, what happens if either of them has no unit, or both have units that are incompatible with each other?

After digging a little bit more, I think the reason why filterwarnings() made me think the indexing tests pass is that when a UnitStrippedWarning is raised (mostly after accessing __array_struct__), the filter replaces it with a DimensionalityError, complaining that it can't convert to dimensionless. This confuses me: when directly accessing an attribute starting with __array_, pytest correctly points out the raised UnitStrippedWarning.

Also, it seems that because DimensionalityError inherits from ValueError, the type of the parameter index of convert_label_indexer is changed from Int64Index to Index (which has dtype="object" ), which does not cast to ndarray, thus keeping the unit information and causing sel() to correctly raise a KeyError. I can reproduce both behaviors using ```python import numpy as np import pytest import xarray as xr

pint = pytest.importorskip("pint") ureg = pint.UnitRegistry()

@pytest.mark.xfail def test_indexing_fail(): array = np.linspace(1, 2, 20) * ureg.m x = np.arange(len(array)) * ureg.s data_array = xr.DataArray(data=array, coords={"x": x}, dims=["x"]) with pytest.raises(KeyError): data_array.sel(x=np.array([13, 5, 9]) * ureg.m)

@pytest.mark.filterwarnings("error:::pint[.*]") def test_indexing_pass(): array = np.linspace(1, 2, 20) * ureg.m x = np.arange(len(array)) * ureg.s data_array = xr.DataArray(data=array, coords={"x": x}, dims=["x"]) with pytest.raises(KeyError): data_array.sel(x=np.array([13, 5, 9]) * ureg.m) ``` Note: the only difference between the two is the decorator.

After digging a bit more into np.prod based on your comment, the situation is more complicated than I thought it was. I forgot that array.prod() has an implementation because pint has a ufunc implementation for it (so, actually np.prod(array) gives a result by way of __array_ufunc__ before https://github.com/hgrecco/pint/pull/764, but it no longer functions that way when __array_function__ is used). I hesitate to delegate to that implementation though, since it gives incorrect units when the axis and/or where arguments are used (see https://github.com/hgrecco/pint/issues/867).

array.any() and array.all() work now because they fall back to .any() and .all() on the magnitude; they aren't explicitly handled (or tested) by pint. Falling back to the magnitude here seems mostly fine, except possibly for offset units. Once a decision is made about the expected behavior in that case, I can add them to the __array_function__ list and add tests in pint.

Based on your tests and examples, I would agree that https://github.com/pydata/xarray/issues/3241 isn't really fixed. I also agree that adding the method tests is a good idea.

Thank you for clarifying about the sel/isel/loc tests. For some reason, I'm unfortunately getting results more consistent with the stripped unit behavior, rather than what your test results are showing. For example:

```python import xarray as xr import numpy as np import pint

unit_registry = pint.UnitRegistry()

array = np.linspace(5, 10, 20).astype(int) * unit_registry.m x = np.arange(len(array)) * unit_registry.s data_array = xr.DataArray(data=array, coords={"x": x}, dims=["x"])

print(data_array.sel(x = [15, 16] * unit_registry.volts)) outputs <xarray.DataArray (x: 2)> <Quantity([8 9], 'meter')> Coordinates: * x (x) int64 15 16 `` whereas it should be giving aKeyError` by your tests, if I'm interpreting them correctly. What happens in your environment when you run the above snippet?

It seems the reason for both the np.median weirdness and the indexing unexpectedly passing is the module level filterwarnings mark. I'm not sure why that happens, but for now I'll remove the module level mark. Also, that means that now all aggregation functions fail for data arrays, so I'll definitely have to add method tests.

np.prod: seems fair, and I don't especially need that function as I have been copying its use from test_sparse.py (though I realize now that there only the method is used). However, using array.prod() works on quantities, so I would expect np.prod to delegate to it. Same with any() and all().

np.median: this also happens in my environment, so you can rule out the environment. Calling it on quantities instead of data arrays works, so the issue is within xarray. Interestingly, outside of pytest: ```python

xr.DataArray(data=np.arange(10).astype(int) * ureg.m).median() <xarray.DataArray ()> <Quantity(4.5, 'meter')> xr.DataArray(data=np.arange(10).astype(int) * ureg.m).median() <xarray.DataArray ()> array(4.5) ``` which seems to be a manifestation of #3241 (which is the reason I now always try different dtypes). I guess that one is not really fixed? Also, would it be better to have a test to check the methods in addition to the numpy functions?

Indexing: when I wrote "indexing works" I meant that my tests for sel(), isel() and loc mostly pass (except from the "calling int / float on scalar results" issue). So when calling sel() and loc with either no units, dimensionless or incompatible units (but existing magnitudes), a KeyError is raised. In the case of both compatible units and existing magnitudes the results are also correct, but scalars trigger the error above. Edit: I have been checking the indexing outside of pytest and there the units are not respected. I have no idea why, but I suspect it is at least related to the problem with np.median. Anyway, what is tested is: ```python

array = np.linspace(5, 10, 20).astype(dtype) * ureg.m x = np.arange(len(array)) * ureg.s data_array = xr.DataArray(data=array, coords={"x": x}, dims=["x"]) with pytest.raises(KeyError): ... data_array.sel(x=[10, 5, 13]) with pytest.raises(KeyError): ... data_array.sel(x=[10, 5, 13] * ureg.dimensionless) with pytest.raises(KeyError): ... data_array.sel(x=[10, 5, 13] * ureg.m)

repr of DataArray fails with DimensionalityError

data_array.sel(x=[10, 5, 13] * ureg.s).data <Quantity([7.63157895 6.31578947 8.42105263], 'meter')> ```

Thank you for the update! Here are responses to each issue brought up:

np.prod (along with np.product and np.nanprod) was not implemented yet since it seems non-trivial to determine the number of unit multiplications that occur given both axis and where arguments in a performant way. I can take a closer look at it though, and hopefully will be able to have an implementation soon. (If anyone has any ideas on this, let me know!)

I passed up np.all and np.any since I associated those with boolean arrays, rather than arrays with units. However, since it sounds like they are needed here, I can add those to the list. (Looking at it though, behavior with non-multiplicative units is unclear right now...see https://github.com/hgrecco/pint/issues/866)

np.allclose should be easy to implement, and I agree it would be nice to not have to have that workaround.

np.maximum and np.minimum appear to be routed through __array_ufunc__ instead of __array_function__. I wasn't aware of that issue until you pointed it out, but looking into the implementation, "maximum" and "minimum" (and any other ufunc not explicitly listed by pint in BaseQuantity.__set_units, BaseQuantity.__copy_units, or BaseQuantity.__prod_units) will be behave in this way since checking input compatibility but not wrapping output is the fallback behavior. This seems like a bug to me that it doesn't raise a UnitStrippedWarning, and at the very least, it is inconsistent with the fail-to-NotImplemented behavior of the __array_function__ implementation. I'll bring up this point at https://github.com/hgrecco/pint/pull/764.

np.argmin/np.argmax/np.sum: there seems to be a lot going on with these at least when I ran the tests locally. First, when called on data_array, they appear to be going deep into xarray/core/duck_array_ops.py and xarray/core/nanops.py rather than ever calling np.nanargmin/np.nanargmax/np.nansum on data_array.data. Second, even if that is the proper pathway, it ends up failing because it ends up calling np.result_type on mixed ndarray and QuantitySequence arguments. This is still unresolved on pint's end (see https://github.com/hgrecco/pint/pull/764#issuecomment-523272038 for the specific issue, https://github.com/hgrecco/pint/issues/845 for the more general discussion). All this also only happens on float dtype...it works just fine with int dtype, which is puzzling to me.

np.median: This one has been driving me crazy. In my environment, when running in pytest, it is erroring out because somewhere np.median being called on data_array is trying to cast the result (a pint QuantityScalar) to a float/int, which is not allowed by pint when the quantity is not dimensionless. But, I have no idea why it is ending up there. However, when replicating the test outside of pytest, it doesn't error, but instead the result is missing the magnitude (with a UnitStrippedWarning raised when going through the __array_struct__ attribute on the quantity. data_array.median(), however, works just fine. I'm not really sure what to do here.

For indexing, when you say "working," would you be able to clarify what your expected behavior is for indexing? Based on https://github.com/pydata/xarray/issues/525#issuecomment-514880353 and the preceding discussion, right now indices will have units stripped, so for me at least, I would expect any attempt at unit-aware indexing to either not work or raise a UnitStrippedWarning.

@jthielen: I tried my tests (those that are currently written) with your PR to hgrecco/pint#764. From the ones I have been using, np.prod, np.all and np.any (and also __matmul__, but you seem to be working on that) are missing. The question is, how many of these should be supported? Also, there might be better ways to compare quantities, but if not: would it be possible to implement np.allclose? I have been using a workaround, but it seems pretty fragile to me. You may be aware of this, but np.minimum and np.maximum and probably others similar to these are working strangely, even with plain pint: if both arguments are quantities, equal units are required, but the result is not a quantity. It correctly fails if the second argument is not a quantity, but works silently if only the second argument is a quantity (still, the result is a plain array).

As I'm not too content with it: could someone have a look at the way they are written? If there are any serious problems, it might be easier to correct them now and I can keep them in mind for new tests.

Lastly, a status report: Assuming there are no bugs in my tests and using andrewgsavage/pint#6, there are some issues with np.argmax/np.argmin, np.median, np.sum (and probably np.prod, but that is currently hidden). Indexing (sel, isel and loc) seems to work except if a single value (e.g. 1 * ureg.m) / index is passed.

ah, okay, that makes sense. I will leave the filterwarning since I believe the warning should not be raised, but mark that test as xfail until a released pint version implements __array_function__.

I noticed you have @pytest.mark.filterwarnings("error::pint.errors.UnitStrippedWarning"), which I believe is what is raising that warning to the error level. What happens when you remove that?

As to why the warning is happening in the first place, I think that is because the __array_function__ implementation in pint right now only exists in a PR. So, when xarray looks for an __array_function__ attribute on the pint Quantity right now, it falls back to this:

https://github.com/hgrecco/pint/blob/2afdc4bf5c5727ed9cef6fdaccb00b88813c4a24/pint/quantity.py#L1438-L1449

Hence, the UnitStrippedWarning. However, hasattr(data, "__array_function__") should still work, albeit for the wrong reason since pint doesn't yet have a true __array_function__ method.

@jthielen: thanks, I will look at these when I have the time.

Could you take a look at the failing test for the DataArray-constructor? It fails because pint emits a UnitStrippedWarning on calling hasattr(data, "__array_function__"), which I believe should not happen (correct me if I'm wrong, though).

@keewis: In case it helps, I've added a bunch of additional __array_function__ implementations in pint in https://github.com/andrewgsavage/pint/pull/6, which will hopefully get merged soon into the main https://github.com/hgrecco/pint/pull/764 PR. If through these tests with xarray you find something missing on pint's end, let me know and I'd be glad to help.

Great! This is going to be exciting. Ping here for any questions!

From what is already done, mypy . fails because of pint (this the same as with quantities in #2956). What would I do to fix that?

V easy - add to setup.cfg like other libraries listed there