I see your point, but xarray could do both, and most accessors I've used come in a pip installable package and we can make that workflow a bit smoother with entry points. IMO it is an advantage of entrypoints that they don't require editing source code, just a 1 line change to a setup.py, setup.cfg, or pyproject.toml.

I wonder how often "users" define their own accessors...I use python functions and modules. The "black magic" you mention breaks most static analysis tooling (type checking, linting, completion) and saves at most a couple of characters, so I never felt the need, but that's a discussion for another day.

@rabernat It seems that great minds think alike ;)

I think the issue is still valid, we just couldn't think of what to name the new API.

@jhamman We have a similar schema package https://github.com/ai2cm/fv3net/tree/master/external/synth, cool to see you confronting the same challenges and advertising your solutions more broadly. One problem we had is that our schema objects ended up being quite verbose: https://github.com/ai2cm/fv3net/blob/master/external/loaders/tests/test__batch/one_step_zarr_schema.json. CDL is a lot more concise.

To be fair, ds.info is not 100% CDL, but it's darn close.

Aren't there xarray extension packages around where this would fit into?

I'm not sure. Any suggestions? Just wondering if xarray has left the door open to this kind of contribution since it 1. already supports other i/o backends 2. creates CDL using ds.info().

I'm sorry too! I don't have any good excuse though...

@jhamman Weren't you talking about an xarray lite (TM) package?

I'm going to close this since I won't be working on it any longer.

NVM, get is already a method. Maybe we could overload it?

Or maybe "get" since it's a synonym of "select" that isn't overloaded with spatial indexing in the code base.

Is there a way in mypy we could use something like overload to specify the above contract here, as an alternative to another method?

That is correct. The output type is predictable from the inputs types. With #4144, mypy may have a chance at detecting this kind of error. Still, I don't know how many users use type-checking. I expect most will only discover this problem at runtime.

It would be better to have an explicit method for subsetting Dataset variables.

I agree. sel_vars is more clear IMO since subset could apply to the coordinates too e.g. a spatial subset.

I’ve run into this as well. It’s not pretty, but my usual work around is to write to a local temporary file and then upload with fsspec. I can never remember exactly which netCDF engine to use...

@crusaderky Thanks for the re-work. For my own benefit, could you explain why that code worked? I remember writing something very similar, and running into mypy errors. My understanding of how mypy intreprets overload seems incomplete.

No problem! I think I am done with this one unless you think its important that I document or test this somehow. Can someone review it?

Okay. Assuming the tests pass, I think this is ready for review. I tried adding a test, but mypy didn't seem to find problems even with code that I know doesn't work (e.g. 'a'+ 1). Is there some strategy for testing tricky type hints like this?

In any case, this code does work: ``` $ cat test_mypy.py (fv3net) import xarray as xr ds = xr.Dataset({"a": ()})

arr = ds['a'] union_obj = ds[['a']]

reveal_locals() $ mypy test_mypy.py (fv3net) test_mypy.py:8: note: Revealed local types are: test_mypy.py:8: note: arr: xarray.core.dataarray.DataArray test_mypy.py:8: note: ds: xarray.core.dataset.Dataset test_mypy.py:8: note: union_obj: Union[xarray.core.dataarray.DataArray, xarray.core.dataset.Dataset] ```

@mathause On further consideration, I think it might not be possible to get this to work. This method has three behaviors: - Mapping -> Dataset - Hashable -> DataArray - else (List): -> Dataset

With my limited understanding of mypy, I think that any two of these is supported by overload, but I'm not sure it's possible to support all 3. I tried several different options, but maybe I am missing something.

Would a good middle ground be something like this? - Hashable -> DataArray - Any -> Union[DataArray, Dataset]

I think this would work since both the input/outputs of the first one are subtypes of the second one. It's not a complete solution, but it would solve the most common problem of ds['a'] returning a union type rather than a DataArray.

~~overload was added in Python 3.5. Is this okay with xarray's backwards compatibility? It would be pretty easy to write a mock overload that does nothing if not.~~

Edit: NVM @mathause's example is in the code already.

@rabernat I learn something new everyday. sorry for cluttering up this PR with my ignorance haha.

@rabernat pointed this PR out to me, and this is great progress towards allowing more database-like CRUD operations on zarr datasets. A similar neat feature would be to read xarray datasets from regions of zarr groups w/o dask arrays.

This variable workaround is awesome @max-sixty. Are there any guidelines on when to use Variable vs DataArray? Some calculations (e.g. fast difference and derivatives/stencil operations) seem cleaner without explicit coordinate labels.

Unfortunately, I don’t think I have much time now to contribute to a general purpose solution leveraging xarray’s built-in indexing. So feel free to add to or close this PR. To be successful, I would need to study xarray’s indexing internals more since I don’t think it is as easily implemented as a routine calling DataArray methods. Some custom numba code I wrote fits in my brain much better, and is general enough for my purposes when wrapped with xr.apply_ufunc. I encourage someone else to pick up where I left off, or we could close this PR.

Thanks so much for the help. This is a good learning experience for me.

That potentially would let you avoid redundant operations on the entire Dataset object.

Yes. This is where I got stuck TBH.

@shoyer Thanks for the comments. I was struggling to incorporate it into Dataset.interp since core.missing is a pretty complicated. Would it be worth refactoring that module to clarify how interp calls are mapped to a given function? Also, most of the methods in interp work like Dataset -> Variables -> Numpy arrays, but the method you proposed above operates on the Dataset level, so it doesn't quite fit into core.missing.interp.

The interpolation code I was working with doesn't regrid the coordinates appropriately, so we would need to do that too.

Ok. I realized this problem occurs only because x was a dimension of the new index idx and the original dataset. Perhaps sel should warn the user or raise an error when this happens.

Ok. I started playing around with this, but I am getting errors when indexing arrays with ND variables. data.sel(x=nd_x) works, but any subsequent operations complain that IndexVariable objects must be 1-dimensional: ```

import xarray as xr import numpy as np npdata = np.tile(np.arange(10), (5, 1)) ... data = xr.DataArray(npdata, dims=['y', 'x'], ... coords={'x': np.r_[:10], 'y': np.r_[:5]}) ... idx = xr.DataArray(npdata, dims=['z', 'x']) ... ... ans = data.sel(x=idx, method='bfill') ... assert set(ans.dims) == {'z', 'y', 'x'} ... print(ans) Traceback (most recent call last): File "<ipython-input-4-1d0ac0cac680>", line 8, in <module> print(ans) File "/Users/noah/workspace/software/xarray/xarray/core/common.py", line 129, in repr return formatting.array_repr(self) File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 463, in array_repr summary.append(repr(arr.coords)) File "/Users/noah/workspace/software/xarray/xarray/core/coordinates.py", line 78, in repr return formatting.coords_repr(self) File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 381, in coords_repr coords, title="Coordinates", summarizer=summarize_coord, col_width=col_width File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 361, in _mapping_repr summary += [summarizer(k, v, col_width) for k, v in mapping.items()] File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 361, in <listcomp> summary += [summarizer(k, v, col_width) for k, v in mapping.items()] File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 307, in summarize_coord coord = var.variable.to_index_variable() File "/Users/noah/workspace/software/xarray/xarray/core/variable.py", line 440, in to_index_variable self.dims, self._data, self._attrs, encoding=self._encoding, fastpath=True File "/Users/noah/workspace/software/xarray/xarray/core/variable.py", line 1943, in init raise ValueError("%s objects must be 1-dimensional" % type(self).name) ValueError: IndexVariable objects must be 1-dimensional ```

So when would interp use this manual interpolation method? Would it try it only if scipy fails or check the dimensions?

After reading the method='ffill' docs, I agree with you.

In my experience, computing w efficiently is the tricky part. The function is slightly different, but metpy uses a lot of tricks to make this work efficiently. A manual for-loop is much cleaner for this kind of stencil calculation IMO. What kind of duck arrays were you thinking of?

I have some numba code which does this for linear interpolation. Does scipy support this pattern?

I am getting the same error too.

Okay. it looks like the docstring for to_unstacked_dataset is in fact incorrect:

https://github.com/pydata/xarray/blob/77a31e56d221245ff7dc10041bf0ca34cab91897/xarray/core/dataarray.py#L1654

I will can submit a PR fixing that soon.

Xarray has a pretty extensive contributor's guide that you might find helpful. In short, the way to contribute changes is to create your own fork of xarray, commit/push some changes, and finally submit a pull request (PR).

Could you link the specific page that is wrong?

Actually, the docs seem correct: https://github.com/pydata/xarray/blob/8da3f67ea583e0588291162067229b2f3ce2993e/xarray/core/dataset.py#L2895

Am I missing something here?

Huh, it looks like that change in the API that happened in the last stages of the PR didn't make it into the example docs. I will try to submit a PR fixing that sometime soon.

Thanks Joe!

phew! Thanks for all the reviews and discussion everyone!

It looks like the CI errors above aren’t related to this PR. There seems to be an issue with to_pandas.

drat. Should be fixed now.

Okay. I responded to @benbovy's last comments and merged the upstream changes to master. How is this looking?

cc @jhamman

is there actual use cases where it is useful to provide multiple dimensions for sample_dims?

Yes! My main use-case does. For example, if you have a weather dataset where every lat/lon pair should be considered a separate "sample", you would use sample_dims=['lat', 'lon'].

It looks like the tests passed. @benbovy How does it look now? Did I fix the issues you mentioned?

It looks like 43834ac8186a851b7 might have fixed this. Let's see if the tests pass.

Does anybody have an idea why the py36-dask-dev tests are failing? None of the errors seems related to this PR.

@benbovy Paradoxically I think making to_stacked_array accept the list of dimensions which won't be stacked is clearer. These unstacked dims needed to shared across all variables, which is a very simple requirement. I realized this after writing a bunch of error-prone validation code. See my last commit for more of my reasoning (5ca9a1d868a). I also updated the docs and tests accordingly.

@benbovy

Is it expected behavior? I guess that not every combination is possible because we're not broadcasting here unlike to_array().

I believe so. In your third case, to_stacked_array basically wants to combine a and b into a DataArray with a dimension y. But that is not possible without broadcasting because b does not have a dimension y. I will try to add a more informative error message.

Thanks for your review @benbovy. I'll try to take a look soon.

Hey @rabernat. Did I resolve your comments?

Would __array_function__ solve the problem with operator precedence? I thought they are separate issues because __mul__ and __rmul__ need not call any numpy functions, and will therefore not necessary dispatch to __array_function__.

To be clear, I think there is some optimal middle ground between the "mega xarray-contrib" package and the current situation. I think the "micro-package" approach works when the collection of micro-packages is being maintained by an active/permanent entity (e.g. Ryan research group). On the other hand, postdocs and grad students are very likely to leave the field entirely within a few years, at which point they will probably stop maintaining their "micro-packages".

I'd also like to thank @teoliphant for weighing in!

Bearing in mind the history of scipy, I agree that the xarray community doesn't need 100% centralization, but there should be some conglomeration. IMO, the current situation of "one graduate student/postdoc per package" is not sustainable.

One additional issue. It seems like pint has some odd behavior with dask. Multiplication (and I assume addition) is not commutative: ``` In [42]: da.ones((10,)) * ureg.m Out[42]: dask.array<mul, shape=(10,), dtype=float64, chunksize=(10,)>

In [43]: ureg.m * da.ones((10,)) Out[43]: dask.array<mul, shape=(10,), dtype=float64, chunksize=(10,)> <Unit('meter')> ```

Hey. I just wanted to bump this PR. How does it look @rabernat @jhamman @shoyer?

@rabernat recent post inspired me to check out this issue. What would this issue entail now that __array_function__ is in numpy? Is there some reason this is more complicated than adding an appropriate __array_function__ to pint's quantity class?

Thanks @rabernat that awesome list looks pretty awesome.

However, I would still advocate for a more centralized approach to this problem. For instance, the NCL has a huge library of contributed functions which they distribute along with the code. By now, I am sure that xarray users have basically reimplemented equivalents to all of these functions, but without a centralized home it is still too difficult to find or contribute new codes.

For instance, I have a useful wrapper to scipy.ndimage that I use all the time, but it seems overkill to release/support a whole package for this one module. I would be much more likely to contribute a PR to a community run repository. I am also much more likely to use such a repo.

I would be more than willing to volunteer for such an effort, but I think it needs to involve multiple people. Various individuals have tried to make such repos on their own, but none seem to have reached critical mass. For example, https://github.com/crusaderky/xarray_extras https://github.com/fujiisoup/xr-scipy I think there should be multiple maintainers, so that if one person drops out, there still appears to be activity on the repo.

Ok. I added an example which is very similar to the DataArray.unstack example.

Okay. It looks like it is passing CI now, so I think it's ready for another look. How does it look? Are you still interested in including this functionality?

@shoyer Sorry. I rebased...

@jhamman I did a little more work on this today. How do you recommend I update this to master? rebase?

Sure, I've been using that as a workaround as well. Unfortunately, that approach throws away all the nice info (e.g. metadata, coordinate) that xarray objects have and requires duplicating much of xarray's indexing logic.

@max-sixty I tend to agree this use case could be outside of the scope of xarray. It sounds like significant progress might require re-implementing core xarray objects in C/Cython. Without more than 10x improvement, I would probably just continue using numpy arrays.

Thanks so much @shoyer. I didn't realize there was that much overhead for a single function call. OTOH, 2x slower than numpy would be way better than 1000x.

After looking at the profiling info more, I tend to agree with your 10x maximum speed-up. A couple of particularly slow functions (e.g. Dataset._validate_indexers) account for about 75% of run time. However, the remaining 25% is split across several other pure python routines.

cc @rabernat

I'd be happy to pick this up again if you think it will go through.

Thanks for the tips everyone. From the dask issue above, it appears that dask.reshape does not respect the boundaries of chunks. Therefore, pulling the first chunk of the stacked array requires pulling several chunks from the original un-reshaped data. This is a bit odd because the first chunk of the stacked array could just be a reshaped version of the first chunk from the unstacked array.

@maxim-lian There is a very short list of such packages hidden in the xarray documention.

In general, there are a ton of these awesome-... repos floating around the internet which just list the useful/related tools/libraries which are related to ... . For example, there are repos out there like awesome-python and awesome-bash. Maybe someone could start an awesome-xarray package.

Sorry for random activity. I accidentally hard reset the master branch nbren12/xarray to pydata/xarray.

cool! Thanks

I just rebased onto master.

no problem. I have always preferred putting operators on new lines though. Didn't realize that was against pep8. oh well 🤷‍♂️ .

Sorry @nbren12 there is still a pep8 error somewhere, and you can freely make it an error as @shoyer suggests.

Hopefully, the commit I just pushed fixes this. I had some global flake8 settings that were messing with my local linting.

Note that currently the test for monotonic coords is for pcolormesh only, while it could be for all 2d plots as well.

I don't think I was having any issues with contour or contourf when I first opened the issue, but it probably does break imshow.

If we start to check input for sanity we might as well raise an error when coordinates are not regularly spaced in the imshow case.

Maybe we could leave this to a later PR.

Yah. I knew my solution before was probably too cute.

On Tue, Feb 6, 2018 at 11:14 AM, Fabien Maussion notifications@github.com wrote:

Thanks! This looks good.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/pull/1885#issuecomment-363532796, or mute the thread https://github.com/notifications/unsubscribe-auth/ABUokvVDiWCnX_JPa5TsDtW8JtQxdTDEks5tSKSZgaJpZM4R4EWU .

There. I think this code should work for all dtypes and 2D coords.

One reason is that it's not obvious if they would like increasing or decreasing coordinates.

For me at least, plt.pcolormesh automatically displays in increasing order even if one of the input arrays is sorted in descending order. This happens all the time with meteorological data available in pressure coordinates (pressure goes down with height). I usually have to manually call plt.gca().invert_yaxis() to flip the y axis. Does xarray.plot's behave differently?

@shoyer That would work with me. Is there any chance people would want to make heatmaps involving categorical variables though?

If we do decide to raise an error, why not go one step further and just sort the coordinates automatically?

On Sun, Feb 4, 2018 at 4:00 PM, Stephan Hoyer notifications@github.com wrote:

@shoyer commented on this pull request.

In xarray/plot/plot.py https://github.com/pydata/xarray/pull/1885#discussion_r165863769:

@@ -750,6 +767,13 @@ def _infer_interval_breaks(coord, axis=0): [ 2.5, 3.5, 4.5]]) """ coord = np.asarray(coord) + + if not _is_monotonic(coord, axis=axis): + warnings.warn("The input coordinate is not sorted in increasing order "

Rather an a warning, why not make this an error? I don't see any use-cases for 2d plots with non-monotonic coordinates. With the current version of xarray, these plots always end up wrong in some way, either by not plotting everything or with bad axis labels.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/pull/1885#pullrequestreview-93863697, or mute the thread https://github.com/notifications/unsubscribe-auth/ABUokuMwNBc8JRgafrF7HNmGfgzxi64Qks5tRkSbgaJpZM4R4EWU .

Ok. I think everything is ready now.

Thanks for your help @fmaussion. Hopefully the commit I just pushed fixes the failing tests. I will work on adding a test and fixing the formatting.

I think automatically sorting is OK for 1D coordinates at least. I agree it is more complicated for other situations.

pcolormesh expects n+1 coordinates

I agree that pcolormesh ultimately works with the mesh corners, but I am pretty sure passing coordinates of the same length also works.

you'll probably have to sort the values beforehand too

True. But we can do it automatically with xarray because it has the coordinate information.

So the question is how much sanity check xarray should do before sending the data to matplotlib, and maybe a warning of some kind would be useful.

I am not sure there there is any circumstance where it would be preferable to plot the scrambled data. Is there some the problem with just adding something like the two lines a wrote above to plot.pcolormesh?

This is pretty easily fixed running python sort_inds = {dim: np.argsort(z[dim].values) for dim in z.dims} z = z.isel(**sort_inds)

I agree that the separate repository model is probably best. However, should it be in just one repository or in many?

Using many repos would solve the domain-specific dependency problem, but the sklearn-contrib packages are not that discoverable IMO. I found two of them via google on separate occasions before realizing that they were part of the same github organization.

Thanks for starting this issue @shoyer. One thing I would be interested to know is how sklearn and tensorflow balance code-quality and API consistency with low barrier to entry. For instance, most of the sklearn contrib packages provide classes which inherit from sklearn's Transformer, BaseEstimator, or Regressor classes, which ensures that all the contrib packages share a common interface.

I would also be very interested in seeing your codes @lamorton. Overall, I think the xarray community could really benefit from some kind of centralized contrib package which has a low barrier to entry for these kinds of functions. So far, I suspect there has been a large amount of code duplication for routine tasks like the fft, since I have also written a function for that.

I guess the main reason barrier for me is initializing the coordinates quickly, which xr.DataArray(np.ones(3, 4), dims=['x', 'y']) doesn't do.

What would use suggest for the signature of a function like xr.ones?

Something like the following would work for me:

```python def ones(shape, dims=None): """Create DataArray of ones with initialized coordinates

Parameters
----------
shape : tuple
    shape of the array
dims : list of str, optional
    list of dimensions with same length as shape. The default will be
    dim_0, dim_1,...,etc
coordinate_initializer : optional
    function which returns the appropriate coordinates. The signature for
    this function must be coordinate_initializer(dim, dim_length). The
    default is to initialize all the coordinates with
    ``np.arange(dim_length)``
"""
pass

```

@shoyer If you are okay with it, I think we might want to leave that to a later date if ever. I am not exactly sure what a useful API for that would be ATM. On the other hand, I have been using the original stack_cat and unstack_cat functions for a couple of months, and they have handled my basic uses pretty well.

For more complicated uses (e.g. taking different subsets of each variable and concatenating the output), I have started working on a project which is similar to sklearn-pandas. Since there are a million ways several xarray variables could be processed/subsetted, stacked and then concatenated, I think this functionality should probably remain in a third party package for now.

Okay. I think I'm done with the updates to the documentation.

I just added some docs, but need to add to whats-new still,

Sorry. I guess I should have made my last comment in the PR.

After using my own version of this code for the past month or so, it has occurred to me that this API probably will not support stacking arrays of with different sizes along shared arrays. For instance, I need to "stack" humidity below an altitude of 10km with temperature between 0 and 16 km. IMO, the easiest way to do this would be to change these methods into top-level functions which can take any dict or iterable of datarrays. We could leave that for a later PR of course.

@rabernat Your point is well taken. I will add some docs/motivation to the reshaping page.

Okay. I just changed the names of the methods and wrote a test case for the problem with the dtype of the stacked dimensions not being preserved by to_stacked_array.

At the moment, I am filling in the missing dimensions with None, so the resulting index has dtype object'. This is then concatenated with the bona fide indices of the variables which are not missing this index, so the index as a whole is cast to the lowest common denominator, which is object. Unfortunately we can't just put Nan in because, NaN is a floating point number thing, so I don't think all numpy dtypes have the equivalent of NaN. @shoyer Do you have any thoughts about how to resolve this?

That naming sounds good to me.

Also, I was having an issue with unstack_cat returning an index with dtype object that I would like to sort out.

Or maybe DataArray.to_dataset_unstacked is good.

Thanks! Yah...I'm not very good at naming things. I think Dataset.to_stacked_array makes sense, but I would typically expect something like from_stacked_array to be a static method of Dataset (e.g. pd.MultiIndex.from_tuples). We could always move unstack_cat to Dataset, but then we miss out on the ability to make sequential method calls.

IMO the behavior of unstack_cat is kind of what I would expect DataArray.to_dataset to do in the case where dim is a stacked coordinate, so maybe we could put it in there.

My understanding is that the "core dimensions" are moved to the end of each input, and then singleton axes are inserted to make them broadcastable in a numpy sense. Is that correct?

Also, it could be clarified what is meant by a core dimension.

I have been playing around with the MultiIndexes as part of #1317, so I could take a stab at implementing some sort of stack which combines the levels etc.

@shoyer I wrote a class that does this a while ago. It is available here: data_matrix.py. It is used like this ```python

D is a dataset

the signature for DataMatrix.init is

DataMatrix(feature_dims, sample_dims, variables)

mat = DataMatrix(['z'], ['x'], ['a', 'b']) y = mat.dataset_to_mat(D) x = mat.mat_to_dataset(y) `` One of the problems I had to handle was with concatenating/stacking DataArrays with different numbers of dimensions---stackandunstackcombined withto_arraycan only handle the case where the desired feature variables all have the same dimensionality. ATM my code stacks the desired dimensions for each variable and then manually callsnp.hstack` to produce the final matrix, but I bet it would be easy to create a pandas Index object which can handle this use case.

Would you be open to a PR along these lines?

Sure. I'd be happy to make a PR once this gets merged. On Sat, Sep 16, 2017 at 10:39 PM Stephan Hoyer notifications@github.com wrote:

Alternatively apply_ufunc could see if the func object has a pre_dask_atop method, and apply it if it does.

This seems like a reasonable option to me. Once we get this merged, want to make a PR?

@jhamman https://github.com/jhamman could you give this a review? I have not included extensive documentation yet, but I am also reluctant to squeeze that into this PR before we make it public API. (Which I'd like to save for another one.)

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/pull/1517#issuecomment-330022743, or mute the thread https://github.com/notifications/unsubscribe-auth/ABUoksuo9P3AJIzemncQQJZ3D5Ga2Opsks5sjLCdgaJpZM4PAViG .