Sorry false alarm the bug that caused salem tests to fail is indeed solved by https://github.com/pydata/xarray/pull/6857

Sorry for the noise.

Is there an older version where it returned copies?

Yes, tests started to fail on salem which I haven't maintained in quite a long time.

Yes sorry I also tested on main and this still occurs. Tagging this as bug.

I've edited a few more sentences, to me this is ready to merge!

I've been struggling with groupby a good deal of my week and I added a warning regarding groupby as well. Feel free to disagree, but I've been unable to get it to work on a large datasets across multiple files yet (see pangeo discourse post).

@Illviljan and other devs - I still think this small change would be a nice addition for accessibility and teaching. Let me know if we can purse or if I should close this.

Side note: with this design (no call to repr() in format()) we would mimick numpy's behavior:

python a = np.array([0.1, 0.2]) a.__format__("") Out[16]: '[0.1 0.2]' a.__repr__() Out[17]: 'array([0.1, 0.2])'

However I'm quite surprised as to why .__format__("") works on non-scalar arrays but not with other specifiers - this is discussed somehow in https://github.com/numpy/numpy/issues/5543 . I think it's okay for xarray to defer to whatever numpy is doing here, but of course we can also keep the status-quo, it's quite a small change after all.

Note that there would be a way to change the behavior to:

python print(da) -> lengthy repr print('{}'.format(da) -> fall back to numpy.__format__

This would break backwards compatibility, but I think it would be my preference - there are some chances that it brakes some code (in documentation pages maybe?), but I don't think that many people rely on '{}'.format(da) to return __repr__ ...

@fmaussion how would you envision this working for non-scalar arrays? Would it just raise?

Yes, similar to numpy:

```python da = xr.DataArray([1, 2, 3]).data print(f'{da}') print(f'{da:d}') [1 2 3]

TypeError Traceback (most recent call last) <ipython-input-93-e6dbd8f1cf00> in <module> 1 da = xr.DataArray([1, 2, 3]).data 2 print(f'{da}') ----> 3 print(f'{da:d}')

TypeError: unsupported format string passed to numpy.ndarray.format ```

The OP bug is now fixed on latest xarray. See #5148 for a related issue with dimension names.

http://xarray.pydata.org/en/stable/generated/xarray.DataArray.rolling.html

What did you try? What did not work out?

I think it's a DeprecationWarning since we're removing it rather than changing its behavior

Ha! I'm no expert either. I learned one thing recently though: DeprecationWarnings are silenced by python per default when not raised directly from __main__ and are therefore invisible for many users.

https://stackoverflow.com/questions/66590950/deprecationwarning-is-not-raised-on-import

https://stackoverflow.com/questions/55377831/difference-between-deprecationwarning-pendingdeprecationwarning-and-futurewarni

Another consequence is that this example in the docs is now meaningless / misleading. I've just come across the series of discussions in cartopy about this issue - it looks quite complex and I'm also not sure yet what to do in xarray in the mean time.

Wouldn't allowing infer_interval=False and passing ordinary arrays as x and y be a good solution?

I'm not sure what you mean by that? You can always set infer_interval=False at the xarray level?

Just adding my +1 here, and also mention that (if memory allows), ds.load() also helps. (related: https://github.com/pydata/xarray/issues/4710)

Thanks for the tip @markelg.

yes, it seems indeed very much related. Closing this in favor of #3961

Thanks a lot for this suggestion! I'll let the others chime in as well, but there are a few things to consider here: - xarray's plotting capability is not meant to be "publication ready" but "as close as possible to the data". I.e. units like degrees_north, etc. have to stay like this because it unambiguously shows where the data came from and how users can access it from the file, etc. - in general, making the plots nicer with opinionated defaults like yours would be a nice addition to the xarray ecosystem. It could live in a separate package, and use the accessor syntax to extend xarray. For example, my own little extension works this way and does opinionated stuff with the accessor: https://salem.readthedocs.io

L67: try to replace

ds_all = xr.merge([xr.open_dataset(f) for f in flist])

by

ds_all = xr.open_mfdataset(flist)

You'll need to install dask for this to work. Note that if you have A LOT of files this might still take some time.

As a not-so-active-but-still-interested xarray dev my opinion doesn't count much, but I would be a proponent of having the rasterio backend live outside of xarray proper.

At the time we wrote the rasterio->DataArray conversion we already noticed a lot of issues regarding differences between the two datamodels, and rio-xarray shows that there is a lot of logic and dev work necessary for this to go smoothly, which would be better handled outside of xarray.

If you want to test builds on RTD before merging you can use the fix-docs branch for that and check its status here: https://readthedocs.org/projects/xray/builds/

Is there a downside to this?

Wouldn't we loose the possibility explore older version of the docs? Or is doctr also providing this service?

It seems so silly to have to reivent readthedocs just because of their CI...

I think we should seriously consider removing all rasterio logic in xarray and make rioxarray an optional dependency to read geotiff files with xarray.

I've noted that the very recent ERA5 data is also cast to float32 at decoding, although the scale and offset params are float64. I think it would be a very valuable addition to xarray to do this properly: @daoudjahdou are you still willing to work on this PR? Otherwise I'll try to pick this up when time permits.

Hi, the docs is here: http://xarray.pydata.org/en/stable/indexing.html#more-advanced-indexing

So your example would read:

```python import numpy as np import xarray as xr da = xr.DataArray(np.arange(16).reshape((4,4)), coords=[('lat', [40,45,50,55]), ('lon', [5,10,15,20])])

points = da.sel(lat=xr.DataArray([40.1, 42.3, 39.78], dims='z'), lon=xr.DataArray([7.1, 6.2, 13.2], dims='z'), method='nearest') ```

Arguably, the "more advanced indexing" title could be renamed to something more explicit.

Can't we have a nightly build of docker images with all the xarray depencencies?

This would help on Travis yes, but not on RTD which doesn't support docker. But otherwise yes, I've also found that CI is faster and more reliable with docker. In terms of resources "in general", I wonder if it wouldn't be nice to share a common testing image base with other packages of the pydata ecosytem.

RTD is struggling with resources, and I can understand that. What is really annoying (and the problem is often the same on travis) is that the conda install is taking a huge part of the build process resources. See e.g. https://github.com/readthedocs/readthedocs.org/issues/6025 , where I ended up using pip for very satisfying results (!).

If I understand well, this is really just about the temporary files, right? In this case, a cleaning block with a :suppress: directive should be enough, yes.

could you change your accessor code to store its state in Dataset.attrs instead?

Yes, although attrs would be bad as well because they are lost in many operations. The current "best practice" would be to use coords for this, as in https://github.com/pydata/xarray/issues/1271#issuecomment-280574680

I never bothered to get this done (my library is quite niche and not mean to become widely used), but if we remove the caching mechanism in xarray, it would be one more incentive to make the switch.

So you work on the assumption that no new potentially useful coords will be added after the first invocation of your accessor?

Yes. In almost all use cases I know of this won't happen. Typically, the users have to create (or open) an xarray object that salem understands before calling the accessor anyway.

Or do you have logic that invalidates your cache every time the state of the coords changes?

No. I didn't think about that before today ;-). But you are right: if users of salem's accessor change the coordinates after calling the accessor, then bad things might happen. Experience shows that this rarely happens (never happened to me), but I can see how this can fire back.

Altogether, these are good arguments to remove the caching mechanism I believe.

Interesting, thanks!

As an accessor maintainer, I can ensure that at least one accessor implementation is storing state ;-). But this state is based on the xarray object itself: for example, we derive georeferencing information and store the names of the coordinate variabless we know are going to be useful to us later. That is, every new call to __init__ based on a modified object will trigger a new parsing, and we don't come into the situation you describe above.

Getting rid of the caching logic would mean some performance loss to us, yes, but I don't know if it's "worse" than the circular reference issue you describe or not.

Yes, network test are not run on travis: https://github.com/pydata/xarray/blob/1a5b6305be4261f04b130a84f7780c9b5dafc943/conftest.py

I'm curious, are there features in rioxarray that could be pushed upstream?

Depends on what the xarray maintainers would like to add. I would definitely like to see the open_rasterio function in rioxarray absorbed back into xarray. Things have just been really slow moving with the other PRs/issues.

Put another way: why don't we put all the logic in rioxarray and make rioxarray an optional dependency of xarray to open rio files?

Thanks @ecarrara !

but the TestRasterio::test_serialization fails

This is not the case anymore, right? This test is here to ensure that the dataarrays we are creating with open_rasterio are compatible with the NetCDF model. I think it is good to try to keep it that way.

This takes about a minute to open for me.

It took me much longer earlier this week when I tried :roll_eyes: Is the bottleneck in the parsing of the coordinates?

Thanks! If you want you can also accept the stackoverflow answer for later reference

But I can't figure out what was the code doing without create this new one "common dimension"

It is doing what is called "orthogonal indexing", but with interpolation. The resulting shape of the output is then (2, 2, 2) in your case, but could be any (t, y, x) as given by the size of each dimension indexer.

Maybe this helps a little: http://xarray.pydata.org/en/stable/indexing.html#vectorized-indexing

Yes, it is possible. It is a bit "less intuitive" at first sight, but powerful and documented here: http://xarray.pydata.org/en/stable/interpolation.html#advanced-interpolation

The call you need to make is:

python blah.interp(longitude=('z', lon), latitude=('z', lat))

How about adding a comment saying that this may/may not work and suggest a workaround

Thanks for the suggestion! I've decided against it because it will be too wordy - let's see if this is really affecting that many people first.

Thanks @grinsted for the PR!

I wonder why it desn't work with pip.

Because of that: https://github.com/mapbox/rasterio#ssl-certs . It will work, but not out-of-the box. I personally don't mind, but we might get a few "your example doesn't work" requests here although it has nothing to do with xarray, this is why I'm asking others devs before merging.

Thanks! Yes, I've tested on a conda env and it works. So my take is that it probably won't work out-of-the box on debian/ubuntu/mint where rasterio is installed via pip.

So my question to @shoyer or @jhamman : should we change this example for the nicer code above, even if a minority of people (pip users like me) might find that this doesn't work for them?

@grinsted what OS are you running this on? I'm a bit concerned about this here: https://github.com/mapbox/rasterio#ssl-certs

I came into this while trying it locally, and even after setting the env variable I still get this error:

RasterioIOError: '/vsicurl/https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif' does not exist in the file system, and is not recognized as a supported dataset name.

I wonder if I'm the only one in that case, or if other users might be affected as well

I'm sorry this has stalled now - I think this is a nice addition.

@jsignell : would you have time to finish this up?

too many references to pandas!

This is a tricky one. In the beginning we wanted people to come from pandas to xarray. It was kind of assumed that everybody uses and loves pandas, therefore will also use and love xarray. (my own experience).

But do we have to build the docs for non-pandas users as well? It could make quite a few things redundant with the pandas documentation...

Thanks! Can you do the same for the other example as well?

(or any other indexing operation)

Correction: only the windowed operations can be done lazily, the other ones will fail unless you load the data first.

Xarray already supports this, the syntax is that of xarray though. It reads the files lazily, then access only what you select, e.g.:

```python import os import urllib.request import matplotlib.pyplot as plt import xarray as xr

Download the file from rasterio's repository

url = 'https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif' urllib.request.urlretrieve(url, 'RGB.byte.tif')

Read the data

da = xr.open_rasterio('RGB.byte.tif') da[:, 400:600, 400:600].plot.imshow(); ```

(or any other indexing operation)

Thanks! @shoyer ! I guess the file I was working with must have been created with nectdf3 then, because otherwise I don't know how they managed to write this out...

It's very python related

I'm saying this because this happened with a real file I was working with. The file had a CLASS attribute on a variable, it can be read by xarray but not written out...

One more request: can you link to a list of disallowed names in the docs somewhere?

I don't know if such a list exists. Do you know? It's very python related, because of setattr. It's confusing that it's only making a problem on variable objects, not on datasets.

So altogether I wouldn't worry much here, this is really an edge case.

This is puzzling. I've edited the example a bit for clarity, but I have no clue what's going on:

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

da = xr.DataArray(np.random.random((3, 4)), dims=["assets", "year"], coords={"year": range(2010, 2018, 2)}, name='var' )

Interpolate directly the DataArray

new_da = da.interp(year=[2011, 2015])

Interpolate via dataset

new_ds = xr.Dataset({"var": da.copy()}).interp(year=[2011, 2015])['var']

Equivalent DAs

xr.testing.assert_identical(new_da, new_ds)

Indexing works

xr.testing.assert_identical(new_da.loc[{"year": 2015}], new_ds.loc[{"year": 2015}])

Assignment doesn't

new_da.loc[{"year": 2015}] = 42 new_ds.loc[{"year": 2015}] = 42 # Fails

IndexError Traceback (most recent call last) <ipython-input-35-65ad561f1692> in <module>() 21 # Assignment doesn't 22 new_da.loc[{"year": 2015}] = 42 ---> 23 new_ds.loc[{"year": 2015}] = 42 # Fails

....

~/Documents/git/xarray/xarray/core/indexing.py in getitem(self, key) 1168 def getitem(self, key): 1169 array, key = self._indexing_array_and_key(key) -> 1170 return array[key] 1171 1172 def setitem(self, key, value):

IndexError: index 2015 is out of bounds for axis 0 with size 2 ```

Done! Didn't edit "what's new" because this is a tiny change.

Should I add tests for this ?

If there is a way, it would be great yes. You can build upon (or adapt) this example here:

https://github.com/pydata/xarray/blob/225868d232219440b188956531d5764ff4cd1b53/xarray/tests/test_backends.py#L3327-L3351

cc @scottyhq who wrote this code, as I don't think any of the xarray team uses this format.

Please also add a line to whats-new.rst

Thanks for the great report!

I can provide this patch in a pull request if needed.

Please, do ;-)

Thanks, I could download them. Can you tell us what the problem with these files is, that we might have to solve in xarray?

Up to now I never thought about that the 'notnull' method is acting on more than only the data itself

All xarray operations will return xarray objects. And xarray will try to match coordinates wherever possible.

However, the coordinates are already mathematically identical

In your example above, they are not. Can you help us to reproduce the error with a Minimal Complete Verifiable Example?

xarray is "coordinate aware", i.e. it will try hard to prevent users doing bad things with non matching coordinates (yes, the fact that your ref and proof are "not entirely consistent somehow regarding coordinates" looks like you are doing bad things ;-).

If I understand what you want, this should do the trick:

python proof["WSS"].where(ref["WSS"].notnull().data) # use .data here to get back to numpy

Thanks!

@benbovy is the checklist still up-to-date? The length of it is a bit scary TBH ;-)

I would be extremely pleased if anyone is willing to jump in and help on the front-end part (HTML/CSS)! See the checklist at the top of this PR. Unfortunately, my limited expertise in this area makes me rather unproductive.

yeah, I guess this is the major issue here. Who could we get in to help out? Does @pydata/xarray know anyone from the extended community with an interest in these things?

A hint in the documentation may hinder others to fall into this pit.

Agreed. Would you like to submit a pull-request?

Also it would be nice to have intervals serialized into netCDF since they are quite common structures.

There are ways to deal with intervals in the CF conventions, but what we really need is a way for xarray to truly understand intervals, which is a much bigger endeavor.

it would save time and frustration if instead of failing, open_dataset would simply issue a warning for any variable it could not decode, and proceed to decode those that it can

I tend to agree, but changing this now would break all code which relies on try ... except logic to check if files can be decoded or not. I'm not sure it's worth changing it at this stage...

@shoyer pep8speaks now support #noqa : https://github.com/OrkoHunter/pep8speaks/issues/96#issuecomment-472310933

I wonder if it's still relevant for us since your change though.

Sorry for being late on this one. Yes, I agree that the "coordinate attribute" is worth a try.

To avoid redundancy in the name one could name the coord variable attributes simply wkt and proj4, like:

<xarray.DataArray 'crs' ()> array(0) Coordinates: crs int64 0 Attributes: wkt: PROJCS["UTM Zone 15, Northern Hemisphere",GEOGCS["WGS 84",D... proj4: +proj=...

I'm going to merge this tomorrow unless there are further concerns

I can't diagnose what's wrong from the error message (something to do with conda it seems)

Some connection error. I restarted Travis, let's see if this happens again.

especially because there is not robust logic as to how to check for coordinates being evenly spaced. I remember this topic has been discussed somewhere already, but I don't remember where.

Thanks @snowman2 ! I don't think there is much discussion necessary about whether we should use WKT or not, the question is rather how? WKT strings are quite complex and verbose, but this should not be a big issue.

We used proj4 strings in xarray because it's what rasterio's to_string would give us back, and we need a canonical string representation for serialization to netCDF. I don't know how many users rely on the rasterio backend's crs attribute to be a proj4 string, but I guess that it's at least a few.

So, what should we do? My preferred solution right now would be to deprecate the crs attribute and replace it with two explicit proj4_crs and wkt_crs attributes to cover all use cases.

Thoughts?

yes, also the current solution isn't too bad as well.

can your merge from master? The failing tests should be solved by https://github.com/pydata/xarray/pull/2720

Maybe we could change it so that attrs is only included if it is non-empty.

How do people feel about this? If folks are fine with it as is, then LGTM.

I'd rather have an empty dict as in the current implementation (like xarray)

It just occurred to me that it would be nice to have some extra info about the data, such as dtype.

II was about to ask about it ;-). What other attribute than dtype were you thinking about?

This is (one more time) an extremely useful addition to xarray - thanks so much @fujiisoup !

I expect some will feel its too soon to merge this

I never wrote py2 code myself so I don't think I'm of the "conservative" type, but I wonder if this PR will make it harder to provide bugfix cherry-picks on the hypothetical py2 compatible 0.11.x branch? Or do we want to ditch this possibility right away and simply focus on the bright py3-only future?

Do we really want a np.masked_array or do we just want an array with np.Nans in it? I prefer nans, but I guess the downside is type conversion for ints?

Masked arrays come with huge performance drops, I don't think we should use these.

I wonder if nodatavals is the only attribute for masking - in the doc (https://rasterio.readthedocs.io/en/latest/topics/masks.html) they also mention nodata, which in my experience is also used

+1 for putting this in "Breaking changes"

Sorry for being so slow. I just done it, will merge tomorrow once the tests pass.

Thanks! Made the change as requested. Regarding the general design:

With this implementation, if one decodes the times, and then saves things back out to a file, a "time bounds" variable will contain units and calendar attributes even though it might not have them initially. I'm not sure if this is a big deal (in this case it is simply adding redundant metadata), but I just wanted to bring it up in case it was a concern.

I agree - I don't think it is a big deal either. It must also be noted that this might break some code in downstream libraries (including mine) which has built workarounds for this and will now error because the bounds are already decoded. Here also, I'm quite confident that this is an edge case but its worth mentioning. Should I put my what's new in "Breaking changes" instead?

I addressed all comments except the fixtures, which seemed a bit overkill (but I simplified the tests). Thanks!

Actually, this shows up in our online docs too

My apologies for the oversight! And thanks for the fix.

Good points!

I do notice the time variable in your file has a bounds attribute, which points to the name time_bnds. Is that something required by CF conventions? We might be able to rely on that.

Yes, bounds is a CF name: http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#cell-boundaries

It seems reasonable to assume that the bounds should share the same units as time, therefore I think that the workaround you are using could actually be implemented in xarray, but the actual implementation might be a bit messy...

(it would be much simpler if CF would prescribe the units to also be available at the the time_bnds level...)

I am actually in favor of 3: also decode time_bnds

Obviously a problem with computation.rst, but could someone help me please?

I doubt this is an xarray problem: we didn't write the documentation with this use case in mind. If you find a way to solve the problem that doesn't brake our Sphinx documentation, then please let us know!

This is the expected behavior, see: http://xarray.pydata.org/en/stable/generated/xarray.DataArray.where.html

but now operates lazily.

If so, how does the dataset access the data if the file is closed? Shouldn't it raise an error instead?

Actually I have to ask: why is it possible to to an operation on a closed dataset? What is happening under the hood which creates this error and wasn't earlier?

Yep, that looks a bug on your end. Were you using autoclose=True before?

No, the exact same piece of code works without error with xarray prev 0.11. Closing this now!

I am slightly regretting not doing a release candidate here,

I should be better placed than anyone to anticipate and test more frequently on master :roll_eyes:

So this reproduces the problem on my machine and travis. Note the misplaced computation out of the with block:

```python import xarray as xr import numpy as np import netCDF4 import os

try: os.remove('test.nc') except OSError: pass

thick = np.zeros((180, 200), np.float32)

ds = xr.Dataset() ds['thick'] = (('y', 'x'), thick) ds.to_netcdf('test.nc')

with xr.open_dataset('test.nc') as ds: dummy = ds.thick.isel(x=('z', [1, 3, 4]), y=('z', [7, 2, 1]))

dummy.min() # this is the issue

with netCDF4.Dataset('test.nc', 'a') as nc: nc.variables['thick'][:] = thick
```

I guess this was actually a bug in our code, but for some reason it worked fine with earlier xarray versions. I'm not sure if this requires action on the xarray side in the end...

I can reproduce this on my machine - also linux. This is going to be hard to track down though.

Did any other library have a warning for this? I think that if we follow the lead of ipython, matplotlib and https://python3statement.org/ then all we have to do is make 0.12 non-installable on py2, and py2 users will stick to 0.11 ad vitam æternam.

It would be quite mean to have them listen to a warning all the time they import xarray on top of that :smile:

ipython docs have a big warning banner though

If not, I can send in a PR to document this bit

I've used the gallery to document these kind of things. It's a bit lost in the doc currently though.

(the use of _mappables is not very obvious; maybe we should make it public?)

Looks like we need it for these kind of feature, right?

Thanks! See https://github.com/pydata/xarray/pull/2498

(Note that GitHub allows to edit files online and does the branching / PR stuff for you, so that small fixes like this can be done exclusively online.)

@am-thyst you are using a relatively old version of python (3.4), and the conda packages are probably not being updated anymore.

Using miniconda is fine, especially if you use conda forge (conda config --add channels conda-forge).

If you want to stay up to date in the future, define a new environment with a newer version of python (3.7 is the latest) and install your packages from there:

conda config --add channels conda-forge conda create -n myenv python=3.7 source activate myenv (linux/mac) or activate myenv (windows) conda install xarray

Doc: https://conda.io/docs/user-guide/tasks/manage-environments.html

Weird that @stickler-ci broke. We're still using it at pandas-gbq and it's working well.

I could try to reactivate it now that I've deactivated it. Still suspicious that clicking @stickler-ci gives a 404 error.

I like pep8speaks for it's less verbose messaging (all in one message that gets updated after each new commit) and for the fact that pandas is also using it, but I have no strong opinion about it.

If we use pep8speaks we should use the same config as pandas though, and we should ask them about #noqa : done here: https://github.com/OrkoHunter/pep8speaks/issues/96

I also have the impression that noqa are ignored by pep8speask, although I haven't seen any related issue on their repo. Probably another reason for a diff_only option for now.

Why do we need diff_only

I was worried about all the unrelated pep8 errors on this PR: https://github.com/pydata/xarray/pull/2463

Although I'm not sure where they come from.

I think we need a set-up file to set diff_only on PR. See the pandas pep8speaks template: https://github.com/pandas-dev/pandas/blob/master/.pep8speaks.yml

Thanks! I just switched off stickler which is still broken anyway.

Appart from a pep8 line length problem LGTM too!

Sitckler is on strike (https://github.com/pydata/xarray/issues/2428)

Closed via https://github.com/rtfd/readthedocs.org/issues/4432

I cannot reproduce this with a dataset with decreasing latitude

Maybe a MPL version issue? Because our docs are also afffected, see: - http://xarray.pydata.org/en/stable/plotting.html#two-dimensions - http://xarray.pydata.org/en/latest/plotting.html#two-dimensions

@fmaussion I'd appreciate it if you could look at this

I will have a look now

Looks like a difficult thing to test, also...

I suppose it would be good to sort out #2422 (Plot2D no longer sorts coordinates before plotting)

That together with https://github.com/pydata/xarray/pull/2398 would be great

But we can save that for another PR :)

See also https://github.com/pydata/xarray/issues/1288 : integrate is the next on my list ;) - I can try to give it a go if @fujiisoup doesn't want to do it himself