Thanks for raising this @vhaasteren ! We want to do what we can to support users from all fields of science :)

I would be okay with that change (especially as it's not really special-casing pint-pulsar, so much as generalizing an existing error-catching mechanism), but would defer to the opinion of @keewis on this.

Solution for those who just found this issue:

Just re-install xarray. pip install -e . is sufficient. Re-installing any way through pip/conda should register the dask chunkmanager entrypoint.

@Illviljan I brought this up in the xarray team call today and we decided that since this only affects people who have previously cloned the xarray repository, are using a development install, and then updated by pulling changes from main; this problem only affects maybe ~10-20 people worldwide, all of whom are developers who are equipped to quickly solve it.

I'm going to add a note into the what's new entry for this version now - if you think we need to do more then let me know.

EDIT: I added a note to whatsnew in https://github.com/pydata/xarray/commit/69445c62953958488a6b35fafd8b9cfd6c0374a5, and updated the release notes.

The only reason I didn't separate the chunkmanager entry points into local and other entry points was simplicity of code.

I didn't realise that might make a difference when it came to whether or not you have to pip install - I assumed that adding a new type of entry point would require re-installing no matter how I implemented it. If that's not the case perhaps we should adjust it (and re-release).

Yes, but I'm wondering what functional difference is that making here?

Have you tried doing the local pip install of the xarray dev version again? I.e. pip install -e . from the xarray folder.

Hmm, it's acting as if dask is not installed/importable. Any idea what's different about your setup vs the xarray CI?

Yes daskmanager is also registered via a different entry point, but that should already be set up to happen by default.

To see which chunk managers it can find you can call

```python from xarray.core.parallelcompat import list_chunkmanagers

list_chunkmanagers() ```

I expect it will return an empty list in your case, but that's the code we should be trying to debug on your system.

np.pad is an interesting example in the context of chunked arrays (xref #6807) - dask implements it (parallelized using various approaches, including map_blocks internally), but cubed currently doesn't implement it because it's not part of the array API standard. (cc @tomwhite)

Closing in favour of #7847

Woooo thanks @dcherian !

@Illviljan thanks for all your comments!

Would you (or @keewis?) be willing to approve this PR now? I would really like to merge this so that I can release a version of xarray that I can use as a dependency for cubed-xarray.

Just merged in https://github.com/pydata/xarray/pull/7820, which should hopefully fix this :crossed_fingers:

If you hover over a node in the SVG representation you'll get a tooltip that shows the call stack and the line number of the top-level user function that invoked the computation. Does that help at all?

That's neat!

When you create a dask graph of xarray operations, the tasks in the graph get useful names according the name of the DataArray they operate on

I realise now that this is not true - but can we make it true for cubed in xarray? Using cubed with xarray creates array's with names like array-002, but couldn't we use the dataarray's .name attribute to give this node of the graph the name "U" for example?

BTW should this be moved to a cubed issue?

I raised it here because it relates to an unfinished part of #7019 - where there is still dask-specific logic for naming individual tasks. I think that to solve this we will need to alter xarray code to allow ChunkManager objects to decide how they want to name their tasks, but using information passed from xarray (i.e. DataArray.name).

(Okay now the failures are from https://github.com/pydata/xarray/pull/7815 which I've separated out, and from https://github.com/pydata/xarray/pull/7561 being recently merged into main which is definitely not my fault :sweat_smile: https://github.com/pydata/xarray/pull/7019/commits/316c63d55f4e2c317b028842f752a40596f16c6d shows that this PR passes the tests by itself.)

I would like to merge this now please! It works, it passes the tests, including mypy.

The main feature not in this PR is using parallel=True with open_mfdataset, which is still coupled to dask.delayed - I made #7811 to track that so I could get this PR merged.

If we merge this I can start properly testing cubed with xarray (in cubed-xarray).

@shoyer @dcherian if one of you could merge this or otherwise tell me anything else you think is still required!

I was not aware of https://github.com/pydata/xarray/issues/6894, which is definitely my bad for not searching properley before setting off smile

No worries! :grin:

It looks like the changes I'm proposing here are probably orthogonal to work in https://github.com/pydata/xarray/issues/6894 though?

I think generally yes they are, I agree.

the goal of this PR is to generalise the existing unit testing to make it a bit easier to run tests with different unit libraries

Any work that helps generalise xarray's support of units beyond specifically just pint is going to be useful!

My main point to draw your attention to is the idea that eventually, one-day, it would be nice to move all array-library specific testing out of the xarray core repo in favour of an approach similar to that proposed in #6894.

I think that testing for unit libraries is a bit less general than the duck array testing stuff, because there's a host of extra information you need to be a unit library compared to a general duck array.

This is also true. Maybe that means for example the base class you are writing here has a long-term future as an optional part of xarray's testing framework in #6894, specifically for use when testing units libraries? Just thinking out loud

Just wanted to drop in and remind people interested in this that we hold a bi-weekly pangeo working group for distributed array computing, which is the perfect place to come and ask about any questions over zoom! I'll be there at 1pm EST today.

Hi @dstansby, thanks for taking initiative on this! Supporting other units-aware packages would be awesome.

Are you aware of our efforts around https://github.com/pydata/xarray/issues/6894? The idea there was to create a general framework for downstream testing of duck-array libraries, including any implementations of units.

I think the ideas you are proposing here are useful and important, but we should probably discuss what we want the end state of duck-array test suites to look like.

cc @keewis

For the benefit of anyone else reading this having come from https://github.com/pydata/xarray/issues/7792 or similar questions - see https://github.com/pydata/xarray/issues/4628 and https://github.com/pydata/xarray/issues/5081 to see what needs to be done. Also see discussion in https://github.com/pydata/xarray/issues/6807 for non-dask lazy backends.

I support this (seems just like what we do for bottleneck) but maybe don't use the word backend for the kwarg again :sweat_smile: In fact as we're only talking about one function could our kwarg literally point to that function? i.e.

python def dot(..., einsum_func=np.einsum): ....

Thanks for raising this @dabhicusp !

So why have that if block at line 396?

Because xarray can wrap many different type of numpy-like arrays, and for some of those types then the self.size * self.dtype.itemsize approach may not return the correct size. Think of a sparse matrix for example - its size in memory is designed to be much smaller than the size of the matrix would suggest. That's why in general we defer to the underlying array itself to tell us its size if it can (i.e. if it has a .nbytes attribute).

But you're not using an unusual type of array, you're just opening a netCDF file as a numpy array, in theory lazily. The memory usage you're seeing is not desired, so something weird must be happening in the .nbytes call. Going deeper into the stack at that point would be helpful.

So while xr.where(cond, x, y) is semantically, "where condition is true, x, else y", da.where(cond, x) is "where condition is true da, else x".

Adding this description to both docstrings would be a helpful clarification IMO.

I'm having problems with ensuring the behaviour of the chunks='auto' option is consistent between .chunk and open_dataset

Update on this rabbit hole: This commit to dask changed the behaviour of dask's auto-chunking logic, such that if I run my little test script test_old_get_chunk.py on dask releases before and after that commit I get different chunking patterns:

python (cubed) tom@tom-XPS-9315:~/Documents/Work/Code/dask$ git checkout 2022.9.2 Previous HEAD position was 7fe622b44 Add docs on running Dask in a standalone Python script (#9513) HEAD is now at 3ef47422b bump version to 2022.9.2 (cubed) tom@tom-XPS-9315:~/Documents/Work/Code/dask$ python ../experimentation/bugs/auto_chunking/test_old_get_chunk.py {'x': (362, 138), 'y': (362, 138)} (cubed) tom@tom-XPS-9315:~/Documents/Work/Code/dask$ git checkout 2022.9.1 Previous HEAD position was 3ef47422b bump version to 2022.9.2 HEAD is now at b944abf68 bump version to 2022.9.1 (cubed) tom@tom-XPS-9315:~/Documents/Work/Code/dask$ python ../experimentation/bugs/auto_chunking/test_old_get_chunk.py {'x': (250, 250), 'y': (250, 250)} (I was absolutely tearing my hair out trying to find this bug, because after the change normalize_chunks became a pure function, but before the change it actually wasn't, so I was trying calling normalize_chunks with the exact same set of input arguments and was still not able to reproduce the bug :angry: )

Anyway what this means is as this PR vendors dask.array.core.normalize_chunks, but the behaviour of dask.array.core.normalize_chunks changed between the version in CI job min-all-deps and the other CI jobs, the single vendored function cannot possibly match both behaviours.

I think one simple way to fix this failure without should be to upgrade the minimum version of dask to >=2022.9.2 (from 2022.1.1 where it currently is).

EDIT: I tried changing the minimum version of dask-core in min-all-deps.yml but the conda solve failed. But also would updating to 2022.9.2 now violate xarray's minimum dependency versions policy?

EDIT2: Another way to fix this should be to un-vendor dask.array.core.normalize_chunks within xarray. We could still achieve the goal of running cubed without dask by making normalize_chunks the responsibility of the chunkmanager instead, as cubed's vendored version of that function is not subject to xarray's minimum dependencies requirement.

I'm having problems with ensuring the behaviour of the chunks='auto' option is consistent between .chunk and open_dataset. These problems appeared since vendoring dask.array.core.normalize_chunks. Right now the only failing tests use chunks='auto' (e.g xarray/tests/test_backends.py::test_chunking_consintency[auto] - yes there's a typo in that test's name), and they fail because xarray decides on different sizes for the automatically-chosen chunks.

What's weird is that all tests pass for me locally but these failures occur on just some of the CI jobs (and which CI jobs is not even consistent apparently???). I have no idea why this would behave differently on only some of the CI jobs, especially after double-checking that array-chunk-size is being correctly determined from the dask config variable within normalize_chunks.

@harshitha1201 if it passes the CI, it should be fine. The error you are getting looks like something to do with importing rasterio, which is an optional backend anyway.

@alrho007 the code for this method on DataArray is in here

https://github.com/pydata/xarray/blob/850156cf80fe8791d45bcaff2da579cffc0cfc35/xarray/core/dataarray.py#L3303

which calls the implementation defined in xarray.core.missing

https://github.com/pydata/xarray/blob/1c81162755457b3f4dc1f551f0321c75ec9daf6c/xarray/core/missing.py#L308

I would start by trying to understand that code (looking at where things are imported to make it work), and then create a small test case example with a monotonically decreasing index that causes a problem. Then try to work out exactly which step in the code causes the issue, and whether it can be generalized to fix the issue.

Thanks @nishtha981 !

I just realised after merging that this PR should in theory have had a corresponding entry in the what's new page, as all PRs are supposed to have.

We won't worry about that this time, but try and remember to add it next time! That way you will also be listed as a contributor on the what's new page.

I've also just told the readthedocs to rebuild just now

Hi @nishtha981 - that's great that you identified what was causing the docs ci builds to fail! I was wondering why that was!

Fixing this for xarray may require more than just re-running the ci. For example it might require us to pin a particular version of a library (here sphinx_book_theme) in order to guarantee the CI works again. If it doesn't work immediately now, what we normally do is to open a new github issue on xarray's issue tracker to track the problem. That way if the problem comes up in multiple PRs we can just refer all of them back to the one issue, until it gets resolved completely.

I've made a final edit to remove the mention of a specific method in the intro, and remove some blanks lines. This looks good to me now so I'll merge it! Thanks @nishtha981 !

In your comment, you said "If you hit commit next to my change, it will merge the change into your pull request " kindly bear with me, I can't seem to find any button that says commit next to your change.

You can ignore this comment - I suggested a change, asked you to commit it, then realised it would be simpler for you if I just used my admin rights to commit it myself, and edited my comment to remove the recommendation. Sorry for the confusion.

Moreover please is there a way I can reach out to you? I tried reaching out to you on Twitter and the provided email, but I have not been successful with that. The mail seems to be invalid.

Sorry! I've been getting a lot of messages from Outreachy applicants and have not managed to reply to them all :cry: The email should not be invalid though - where did you send it?

If you have questions about xarray itself they should be raised on the repository though - you don't need to contact me directly for that. Also if you raise them publicly on the repository it gives other people a chance to answer if I don't see it.

If you have questions about Outreachy specifically then ask them either on the discord channel or by emailing me (thomas dot nicholas at columbia dot edu).

Hi @harshitha1201 - thanks for this!

We do already have a section covering these methods in https://docs.xarray.dev/en/stable/user-guide/computation.html#missing-values. I suggest that we don't need to duplicate all of this information on the FAQ page.

That said the examples and explanation you have written here are still useful! Perhaps they can either be used to improve the page I just linked, or go into the docstrings of those particular methods.

For the FAQ page instead I think we probably just want to provide a summary in a couple of sentences and a link to the more detailed information on specific methods. The summary should mention that - xarray can handle missing values, - it uses np.NaN to do so, - most computation methods will automatically handle missing values appropriately, - aggregation methods have a skipna argument, - plotting will just leave them as blank spaces (link to plotting page), - we have a set of special methods for manipulating missing and filling values (link to here).

Good idea @dcherian

Could you specifiy what kind of function you would like the bot to perform?

When a github user creates their first ever issue or pull request on the xarray repository, the bot would reply to them welcoming them with a friendly message, thanking them for their interest, and point them towards useful links like contributing guidelines and contact channels.

:+1: to not having to write to_dataset before saving things.

I think the .to_X vs e.g. .save(engine='X') is important enough that we should think about it a bit more. The response to #7496 doesn't seem conclusive to me (yet).

A twitter poll? Then a discussion in a meeting?

Hi @Amisha2778 - unfortunately this issue has already been resolved by pull request #7625. We just forgot to close this issue once it was resolved, sorry about that! (The closing normally happens automatically but apparently didn't this time.)

You are welcome to work on any other issue you like, but I will close this one now.

Hi @Amisha2778 - great to hear you are interested. You don't need my permission - please have a go at solving any issue that looks interesting to you, and please ask questions if you have any difficulties!

Thanks @harshitha1201 ! A great usability improvement.

I've made a bare-bones cubed-xarray package to store the CubedManager class, as well as any integration tests (yet to be written). @tomwhite you should have an invitation to be an owner of that repo. It uses the entrypoint exposed in this PR to hook in, and seems to work for me locally :grin:

Thanks for the review @dcherian! I agree with basically everything you wrote.

The main difficulty I have at this point is non-reproducible failures as described here

I would like to get to the point where you can use xarray with a chunked array without ever importing dask. I think this PR gets very close, but that would be tricky to test because cubed depends on dask (so I can't just run the test suite without dask in the environment

I just released Cubed 0.6.0 which doesn't have a dependency on Dask, so this should be possible now.

Actually testing cubed with xarray in an environment without dask is currently blocked by rechunker's explicitly dependency on dask, see https://github.com/pangeo-data/rechunker/issues/139

EDIT: We can hack around this by pip installing cubed, then pip uninstalling dask as mentioned here

Does this mean my comment https://github.com/pydata/xarray/pull/7019#discussion_r970713341 is valid again?

Yes I think it does @headtr1ck - thanks for the reminder about that.

I now want to finish this PR by exposing the "chunk manager" interface as a new entrypoint, copying the pattern used for xarray's backends. That would allow me to move the cubed-specific CubedManager code into a separate repository, have the choice of chunkmanager default to whatever is installed, but ask explicitly what to do if multiple chunkmanagers are installed. That should address your comment @headtr1ck.

This came up in the xarray office hours today, and I'm confused why this PR made any difference to the behavior at all? The .data property just points to ._data, so why would it matter which one we check?

Thanks @dcherian ! Once I copied that explicit indexer business I was able to get serialization to and from zarr working with cubed!

```python In [1]: import xarray as xr

In [2]: from cubed import Spec

In [3]: ds = xr.open_dataset( ...: 'airtemps.zarr', ...: chunks={}, ...: from_array_kwargs={ ...: 'manager': 'cubed', ...: 'spec': Spec(work_dir="tmp", max_mem=20e6), ...: } ...: ) /home/tom/Documents/Work/Code/xarray/xarray/backends/plugins.py:139: RuntimeWarning: 'netcdf4' fails while guessing warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning) /home/tom/Documents/Work/Code/xarray/xarray/backends/plugins.py:139: RuntimeWarning: 'scipy' fails while guessing warnings.warn(f"{engine!r} fails while guessing", RuntimeWarning)

In [4]: ds['air'] Out[4]: <xarray.DataArray 'air' (time: 2920, lat: 25, lon: 53)> cubed.Array<array-004, shape=(2920, 25, 53), dtype=float32, chunks=((730, 730, 730, 730), (13, 12), (27, 26))> Coordinates: * lat (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0 * lon (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0 * time (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00 Attributes: GRIB_id: 11 ...

In [5]: ds.isel(time=slice(100, 300)).to_zarr("cubed_subset.zarr") /home/tom/Documents/Work/Code/xarray/xarray/core/dataset.py:2118: SerializationWarning: saving variable None with floating point data as an integer dtype without any _FillValue to use for NaNs return to_zarr( # type: ignore Out[5]: <xarray.backends.zarr.ZarrStore at 0x7f34953033c0> ```

@Ravenin7 do you need some guidance on how to add a test for this? Happy to help if so. It would be great to get this fix merged soon because a number of other people have encountered the same bug!

I tried opening a zarr store into xarray with chunking via cubed, but I got an error inside the indexing adapter classes. Somehow the type is completely wrong - would be good to type hint this part of the code, because this happens despite mypy passing now.

```python

create example zarr store

orig = xr.tutorial.open_dataset("air_temperature") orig.to_zarr('air2.zarr')

open it as a cubed array

ds = xr.open_dataset('air2.zarr', engine='zarr', chunks={}, from_array_kwargs={'manager': 'cubed'})

fails at this point

ds.load() ```

This still works fine for dask.

Thanks @tomwhite - I think it might make sense for me to remove the CubedManager class from this PR and instead put that & cubed+xarray tests into another repo. That keeps xarray's changes minimal, doesn't require putting cubed in any xarray CI envs, and hopefully allows us to merge the ChunkManager changes here earlier.

Places dask is still explicitly imported in xarray

There are a few remaining places where I haven't generalised to remove specific import dask calls either because it won't be imported at runtime unless you ask for it, cubed doesn't implement the equivalent function, that function isn't in the array API standard, or because I'm not sure if the dask concept used generalises to other parallel frameworks.

• [ ] open_mfdataset(..., parallel=True) - there is no cubed.delayed to wrap the open_dataset calls in,
• [ ] Dataset.__dask_graph__ and all the other similar dask magic methods
• [ ] dask_array_ops.rolling - uses functions from dask.array.overlap,
• [ ] dask_array_ops.least_squares - uses dask.array.apply_along_axis and dask.array.linalg.lstsq,
• [ ] dask_array_ops.push - uses dask.array.reductions.cumreduction

I would like to get to the point where you can use xarray with a chunked array without ever importing dask. I think this PR gets very close, but that would be tricky to test because cubed depends on dask (so I can't just run the test suite without dask in the environment), and there are not yet any other parallel chunk-aware frameworks I know of (ramba and arkouda don't have a chunks attribute so wouldn't require this PR).

@alrho007 great! We would welcome a pull request to update this :)

Thanks @remigathoni ! We would welcome a pull request that made these improvements :slightly_smiling_face:

Hi @mahamtariq58, thanks for your interest! We don't normally assign issues to individuals, you are just welcome to have a go at solving any issue that interests you.

Hi @alrho007 - thanks for your interest in contributing to xarray!

The code for the string accessor is in xarray/core/accessor_str.py. That's where you would need to make changes.

The idea would be to wrap xr.DataArray.str.mod in xr.DataArray.str.format.

@mathause could you expand on what you mean exactly? I'm new to this part of the codebase. There is already a .format method, are you talking about changing its behaviour?

xr.DataArray([0, 1, 2]).str would raise an error right since it's int type?

Also @ahuang11 this seems fine? Because the .str accessor doesn't actually perform a check that the passed data is a str type immediately?

Oh right my bad - I needed to click the lower of the two "View Docs" buttons :sweat_smile:

That's definitely confusing, so I'll add this screenshot to this PR.

Wait hang on the docs build in the CI is not pointing to my PR - it's pointing to the stable branch! That's not right.

I'm making progress with this PR, and now that @tomwhite implemented cubed.apply_gufunc I've re-routed xarray.apply_ufunc to use whatever version of apply_gufunc is defined by the chosen ChunkManager. This means many basic operations should now just work:

```python In [1]: import xarray as xr

In [2]: da = xr.DataArray([1, 2, 3], dims='x')

In [3]: da_chunked = da.chunk(from_array_kwargs={'manager': 'cubed'})

In [4]: da_chunked Out[4]: <xarray.DataArray (x: 3)> cubed.Array<array-003, shape=(3,), dtype=int64, chunks=((3,),)> Dimensions without coordinates: x

In [5]: da_chunked.mean() Out[5]: <xarray.DataArray ()> cubed.Array<array-006, shape=(), dtype=int64, chunks=()>

In [6]: da_chunked.mean().compute() [cubed.Array<array-009, shape=(), dtype=int64, chunks=()>] Out[6]: <xarray.DataArray ()> array(2) ```

(You need to install both cubed>0.5.0 and the main branch of rechunker for this to work.)

I still have a fair bit more to do on this PR (see checklist at top), but for testing should I:

1. Start making a test_cubed.py file in xarray as part of this PR with bespoke tests,
2. Put bespoke tests for xarray wrapping cubed somewhere else (e.g. the cubed repo or a new cubed-xarray repo),
3. Merge this PR without cubed-specific tests and concentrate on finishing the general duck-array testing framework in #6908 so we can implement (b) in the way we actually eventually want things to work for 3rd-party duck array libraries?

I would prefer not to have this PR grow to be thousands of lines by including tests in it, but also waiting for #6908 might take a while because that's also a fairly ambitious PR.

The fact that the tests are currently green for this PR (ignoring some mypy stuff) is evidence that the decoupling of dask from xarray is working so far.

(I have already added some tests for the ability to register custom ChunkManagers though.)

There is a link in the docs to an xarray wiki:

If your test requires working with files or network connectivity, there is more information on the `testing page <https://github.com/pydata/xarray/wiki/Testing>`_ of the wiki.

This link is broken, but I'm not sure what to replace it with. I didn't even know xarray had a wiki!! Does it still exist?

Thanks for reporting @arfriedman .

The link is supposed to be referring to the section immediately following the Datetime Indexing subsection. I think the underscores in datetime_component_indexing might be what's breaking the link.

If anyone wants to make their first small contribution to xarray this would be a great start!

Wow thanks for reporting this @Metamess ! Do you see a simple fix within merge_core?

See #3980 - this issue should be solved as part of #3980, but see that issue for more general discussion of subclassing.

Closing this as having answered the original question. If anyone wants to discuss mosaicing rasters in more detail we can raise another issue for that.

Or just xr.concat.

The chunk_by_labels functionality seems quite useful even when not talking about times, so I would be :+1: for that kind of option.

On the API question is there anywhere else in xarray where we have made some choice about how to let the user choose between specifying via indexes or labels? Apart from just .isel vs .sel I mean

Or indeed at least make it clearer in the docs that something like drop_indexes or reset_coords should be used first in order to skip auto-alignment for some variables.

I think we should do this regardless. I don't know of anywhere in the docs that these kind of subtleties with concat are clearly documented.

I guess easiest for a concat version with no auto-alignment would be to drop the index when such case happens.

Right - in this case that would have given the intuitive result.

We could get halfway to a better xr.concat by changing https://github.com/pydata/xarray/issues/2064 IMO.

I propose join="exact", data_vars="minimal", coords="minimal", compat="override"

That wouldn't have helped with this specific issue though right?

it was triggering a load

Can we not just test the in-memory performance by .load()-ing first? Then worry about dask performance? That's what I was vaguely getting at in my comment, trying the in-memory performance but also plotting the dask graph.

Could you show the example that's this slow, @TomNicholas ? So I can play around with it too.

I think I just timed the difference in the (unweighted) "real" example I gave in the notebook. (Not the weighted one because that didn't give the right answer with flox for some reason.)

One thing I noticed in your notebook is that you haven't used chunks={} on the open_dataset. Which seems to trigger data loading on strange places in xarray (places that calls self.data), but I'm not sure this is your actual problem.

Fair point, worth trying.

We'll fix the compatibility issues, but we first need to understand what the expectations on something like data.shape should be in these circumstances.

At a minimum, xarray expects .shape, .ndim and .dtype to always be defined. (And the number of dims to match the shape, which Joe's example above implies aesara doesn't do?) On top of that there are extra expectations about slicing and broadcasting changing shape in the same ways as it does for numpy arrays. (@keewis correct me if I've mis-stated this or missed something important here!)

For a shared variable, it's always possible to get a value for data.shape by referencing the underlying data, but the reason we don't do that by default is—in part—due to the fact that shared variables can be updated with values that have different shapes (but the same dtypes and number of dimensions).

This sounds a bit similar to discussions we have been having about wrapping ragged arrays in xarray, for which there are multiple ways you might choose to define the shape.

The simplest way to guarantee that aesara can be wrapped by xarray is if aesara conformed to the array API standard, and they have a test suite you can use to check that. We are also working on our own testing framework that duck-typed array libraries like aesara could import to quickly test integration with xarray.

Q: Use xhistogram approach or flox-powered approach?

@dcherian recently showed how his flox package can perform histograms as groupby-like reductions. This begs the question of which approach would be better to use in a histogram function in xarray.

(This is related to but better than what we had tried previously with xarray groupby and numpy_groupies.)

Here's a WIP notebook comparing the two approaches.

Both approaches can feasibly do: - Histograms which leave some dimensions excluded (broadcast over), - Multi-dimensional histograms (e.g. binning two different variables into one 2D bin), - Normalized histograms (return PDFs instead of counts), - Weighted histograms, - Multi-dimensional bins (as @aaronspring asks for above - but it requires work - see how to do it flox, and my stalled PR to xhistogram).

Pros of using flox-powered reductions:

• Much less code - the flox approach is basically one call to flox.
• Fewer codepaths, with groupby logic and all histogram functionality flowing through the flox.xarray_reduce codepath.
• Likely clearer code than the kinda impenetrable reshaped bincount logic lurking in the depths of xhistogram.
• Supporting new features (e.g. multidimensional bins) should be simpler in flox because the options don't have to be propagated all the way down to the level of the np.bincount caller.

Pros of using xhistogram's blockwise bincount approach:

• Absolute speed of xhistogram appears to be 3-4x higher, and that's using numpy_groupies in flox. Possibly flox could be faster if using numba but not sure yet.
• Dask graphs simplicity. Xhistogram literally uses blockwise, whereas the flox graphs IIUC are blockwise-like but actually a specially-constructed HLG right now. (Also important for supporting other parallel backends.) I suspect that in practice both perform similarly well after graph optimization but I have not tested this at scale, and flox's graph might be more sensitive to extra steps in the calculation like adding weights or normalising the result.

Other thoughts:

• Flox has various clever schemes for making general chunked groupby operations run more efficiently, but I don't think histogramming would really benefit from those unless there is a strong pattern to which values likely fall in which bins, that is known a priori.
• Deepak's example using flox uses pandas.IntervalIndex to represent the bins on the result object, whereas xhistogram just returns the mid-points of the bins, throwing that info away. This seems like a cool idea on it's own, but probably requires some extra work to make sure it's handled by the indexes refactor and the plotting code.
• In my comparison notebook here's something I'm missing that's causing my "real example" (from xhistogram docs) to not actually use the provided weights. I suspect its something simple, any idea @dcherian?

xref https://github.com/xgcm/xhistogram/issues/60, https://github.com/xgcm/xhistogram/issues/28

(I was already typing this when Ryan posted so I'll finish anyway :sweat_smile:)

To clarify, what @jhamman is suggesting in this specific issue is xarray wrapping aesara, as opposed to aesara wrapping xarray.

Both of these goals would be great, but in xarray we would particularly love to be able to wrap aesara because: - We already have a large community of users who use xarray objects as their top-level interface, - It continues a larger project of xarray being generalized to wrap any numpy-like "duck array". - Our ultimate goal is for xarray users to be able to seamless switch out their computing backend and find which library gives them the best performance without changing the rest of their high-level code.

For xarray to wrap aesara aesara needs to provide a numpy-like API, ideally conforming to the python array api standard. If aesara already does this then we should try out the wrapping right now!

If you're interested in this topic I invite you to drop in to a meeting of the Pangeo working group on distributed arrays! We have so far had talks from distributed computing libraries including Arkouda, Ramba, and cubed, all of which we are hoping to support as compute backends.

If anyone is also separately interested in using xarray inside PyTensor / aesara then that's awesome, but we should try to track efforts in that direction on a different issue to keep this distinction clear. We plan to better support that direction of wrapping soon by fully exposing our (semi-private internal currently) lightweight Variable class.

Thanks @dcherian!

Wait what version are we on right now?? :sweat_smile: Are you in middle of issuing a release? The whats-new in main doesn't match what's actually been released atm

I think this PR is ready now, it just fails mypy (@Illviljan I added a py.typed file to datatree but the xarray mypy CI is still not happy with it).

Once this is merged we can push on with implementing open_datatree (https://github.com/pydata/xarray/pull/7437 @jthielen) and I can take a crack at fixing https://github.com/xarray-contrib/datatree/issues/146 in xarray in a follow-up PR.

@khider we are more than happy to help with digging into the codebase! A reasonable place to start would be just trying the operation you want to perform, and looking through the code for the functions any errors get thrown from.

You are also welcome to join our bi-weekly community meetings (there is one tomorrow morning!) or the office hours we run.

Hi @khider , thanks for raising this.

For those of us who haven't tried to use non-nanosecond datetimes before (e.g. me), could you possibly expand a bit more on

However, most of the interesting functionalities of xarray don't seem to support this datetime out-of-box:

specifically, where are errors being thrown from within xarray? And what functions are you referring to as examples?

@pydata/xarray should this feature be added to our development roadmap? It's arguably another approach to making more flexible data structures...

@drtodd13 mentioned today that ramba doesn't actually require explicit chunks to work, which I hadn't realised. So forcing wrapped libraries to implement an explicit chunks method might be too restrictive. Ramba could possibly work entirely through the numpy array API standard.

@drtodd13 tagging you here and linking my notes from today's distributed arrays working group meeting for the links and references to this PR.

This seems like a good idea to me, but I don't know much about this part of the codebase. We should at minimum state in the docstring of these method whether they are required or optional.

From the errors thrown in the tests it seems set_variable is not defined by backends quite regularly at least.

But then da.data will be of this protocol type and not the array class that you assume it has. For internal type checking this is what we want but for the user this will be confusing.

When will the user be using this type annotation? Isn't all this typing stuff basically a dev feature (internally and downstream)?

Doesn't the python array api standard effort have some type of duck array protocol we could import? I feel like this has been mentioned before. Then we would start with @Illviljan 's suggestion and replace it with the correct duck array protocol later.

We might also consider that in the context of different distributed array backends dask arrays define a more specific API that includes methods like .chunk() too. Standardisation of that is a long-term dream though, not an immediate problem.

Thanks for these comment @paigem! It's useful to hear things that might not occur to us.

Good points @Illviljan. Perhaps we should restructure to make it clear that running tests / pre-commit locally is optional, but still show people how to do it?

@ronygolderku thanks for your example. Looks like it fails for the same reason as was mentioned for some of the other examples above.

Holla when you want a review :)

Thanks for the suggestion @nbren12 !

Whilst I agree that this would be a more "correct" way of providing accessor functionality, I think there is a big downside in that Entrypoints are quite a lot harder for users to use than the accessors are.

Whilst the way the accessors actually add the method is kind of black magic, all the user has to do is copy-paste the example from the docs and change the methods and accessor name to what they want, and it will immediately work.

Adding an entrypoint requires going into your setup.py (which you wouldn't even have if you're running a script or notebook), and it's conceptually complicated. The accessors stuff is useful even for relatively novice users, so I don't really want to make it harder to understand...

Interested in what other people think though?

(I'm looking for issues to tag as beginner-friendly)

Shall we also add a top-level import: xr.show_backends like xr.show_versions. This would be a lot more discoverable.

This is easy right? You could just add an xr.show_backends that points to the same function as xr.backends.list_engines? And point out in the docstring that they are actually the same function?

@gabicca given that we do say in our docs that we don't recommend subclassing xarray objects, I'm going to close this issue as a downstream "subclass at your own risk" problem. If you would like it to be easier to subclass, we would welcome your input here.

The datatree class you mentioned is very interesting, is it yet available for me to use as a package?

Also yes you can find datatree here, though it's a bit experimental compared to the rest of xarray.

@NilsHoebe to clarify, are you still unable to open the desired group even when you use group as the kwarg? If so a minimal example would help us to help you!

Integrating upstream into xarray might also help with people trying to open their nested data formats as a datatree objects, because then we can immediately begin integrating with xarray's backend engines.

See for example this datatree issue asking about opening grib files as a datatree. It would be nice to be able to do

open_datatree("data.grib", engine="cfgrib")

Would it mean that if someone wants to later add any feature "x" or "y" into Xarray, they just need implementing the feature for Dataset (and possibly DataArray) and it will be guaranteed to work with Datatree?

Basically yes, it would immediately work with Datatree. Datatree currently implements most dataset methods by literally copying them and their docstrings, and they work by mapping the method over every node in the tree. We could integrate Datatree in such a way that the additional developer effort to get a method on dataset working on Datatree would be negligible (think adding a single element with the method name to an internal list, or copy-pasting a docstring).

I don't think it is ideal to have such non-synchronized state within Xarray itself.

This is an argument for waiting before integrating.

I'm just speaking generally from my experience of having struggled while doing some heavy refactoring in Xarray recently :)

I appreciate the input @benbovy! I think the main difference between this effort and your (heroic) indexes effort is that Datatree doesn't touch any existing API.

I guess my main concern is that integrating prematurely into Xarray might give a false sense of stability - I don't want to later realize I should redesign Datatree, and have people be annoyed because they thought it was as stable as the rest of xarray.

Why? Are its dependencies different from Xarray?

No, datatree has no additional dependencies. I was just asking because if we went for the "import from second repository" plan we may want to test that the import works as part of our CI. Not a major issue though.

There is also at least one bug in datatree that cannot be fixed without a (small) change to xarray, and having datatree as an optional import means I could fix it here.

Thanks for this bug report. FWIW I have also seen this bug recently when helping out a student.

The question here is whether this is an xarray, numpy, or a netcdf bug (or some combo). Can you reproduce the problem using to_zarr()? If so that would rule out netcdf as the culprit.

That's a useful observation, thank you @maawoo!

This comes from the way we generate our code for the many different aggregations xarray can perform. We actually use this script to automatically generate all the source code for all the aggregations in this file. That script has a template that is filled in for each method.

Currently the template looks like this

python TEMPLATE_REDUCTION_SIGNATURE = ''' def {method}( self, dim: Dims = None, *,{extra_kwargs} keep_attrs: bool | None = None, **kwargs: Any, ) -> {obj}: """ Reduce this {obj}'s data by applying ``{method}`` along some dimension(s). Parameters ----------''' where in the case of variance the method is just var so "variance" isn't in the generated docstring anywhere.

How might we fix this? One immediate thought that might help is to change the template to use a method_name and a long_name, where method_name is var but long_name is variance for example. This shouldn't be particularly difficult, and we would welcome a PR if you would be interested in contributing? We would help you out :slightly_smiling_face:

Or we might change the docstrings in some other, more granular way. Adding examples to aggregation methods would also have to deal with the fact they are autogenerated https://github.com/pydata/xarray/issues/6793

Using xr.Grouper has the advantage that you don't have to start guessing about whether or not the user wanted some complicated behaviour (especially if their input is slightly wrong somehow and you have to raise an informative error). Simple defaults would get left as is and complex use cases can be explicit and opt-in.

I kinda think correctness by default is more important than performance, especially if the default performance isn't prohibitively slow.

@keewis https://github.com/pydata/xarray/pull/7338 mentions datatree, so we could refer people back to there from now on?

I realise that this was probably a function intended for internal use only

Yes, everything under xarray.core is internal, and unfortunately internal functions aren't publicly supported, so we provide no guarantee that they will still work/exist between versions.

Is there a better way to do this? What replaces this function?

@benbovy should know!

The email that you have listed here doesn't work (bounced back).

Oops - use thomas dot nicholas at columbia dot edu please!

I should be able to join today as well @jpivarski ! Will need the zoom address

Hi @airton-neto - here is a much shorter example that reproduces the same error

```python url = 'https://rda.ucar.edu/thredds/dodsC/files/g/ds084.1/2022/20220201/gfs.0p25.2022020100.f000.grib2'

dataset = xr.open_dataset(url, cache=True, engine="netcdf4") python from dagster._utils import frozenlist

variables = list(['u-component_of_wind_height_above_ground']) # <--- This way it runs

variables = frozenlist(['u-component_of_wind_height_above_ground'])

dataset[variables] ```

```python

KeyError Traceback (most recent call last) File ~/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/core/dataset.py:1317, in Dataset._construct_dataarray(self, name) 1316 try: -> 1317 variable = self._variables[name] 1318 except KeyError:

KeyError: ['u-component_of_wind_height_above_ground']

During handling of the above exception, another exception occurred:

KeyError Traceback (most recent call last) Input In [43], in <cell line: 1>() ----> 1 dataset[variables]

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/core/dataset.py:1410, in Dataset.getitem(self, key) 1408 return self.isel(**key) 1409 if utils.hashable(key): -> 1410 return self._construct_dataarray(key) 1411 if utils.iterable_of_hashable(key): 1412 return self._copy_listed(key)

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/core/dataset.py:1319, in Dataset.construct_dataarray(self, name) 1317 variable = self._variables[name] 1318 except KeyError: -> 1319 , name, variable = _get_virtual_variable(self._variables, name, self.dims) 1321 needed_dims = set(variable.dims) 1323 coords: dict[Hashable, Variable] = {}

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/xarray/core/dataset.py:171, in _get_virtual_variable(variables, key, dim_sizes) 168 return key, key, variable 170 if not isinstance(key, str): --> 171 raise KeyError(key) 173 split_key = key.split(".", 1) 174 if len(split_key) != 2:

KeyError: ['u-component_of_wind_height_above_ground'] ```

I'm not immediately sure why the elements of the list are not properly extracted, but the frozenlist class is a dagster private internal, and the problem is specific to how that object interacts with the current xarray code. As a downstream user of xarray, can dagster not just change it's code to pass the type xarray expects (i.e. a normal list)?

Having said that if you want to submit a PR which fixes this bug (and doesn't require special-casing dagster somehow) then that would be welcome too!

Well that's frustrating - I didn't see that.

They also mention that if PEP 637 hadn't been rejected there would be more use, which is ironic because we also wanted PEP 637 for indexing (like da[time=slice(5, 10)]).

Big moves!!

If one wants to get fancy, we could wrap the list_engines dict into a custom dict with a nice repr (maybe even html repr).

Whilst I don't want to discourage eager devs, in one of the community calls we did discuss how this would almost certainly be overkill for a little-used advanced feature.

LGTM as long as it renders nicely.

It renders nicely locally, but any idea why the RTD build fails so opaquely?