Thanks @benbovy - with explicit methods now to produce the result with or without index, I think we can close this now :smiley:

I didn't find any solution in xarray - I think I ended up just dropping the conflicting variable...

The vote seems to be for a separate cumulative_integrate method - I've made that change.

One thing that could be nice is an example in http://xarray.pydata.org/en/stable/plotting.html - up to you.

It seemed odd to add a surface example when contour and countourf weren't shown - so I've added a new subsection for 'Other types of plot' and shown all three.

Thanks @mathause - master merged now.

Thanks @mathause! I've used your suggestion to turn my MVCE into a test which fails on master, but passes with the fix in this PR.

There's a cumtrapz in https://github.com/fujiisoup/xr-scipy/blob/master/xrscipy/integrate.py. Does that help?

I did a quick check and it seems like applying scipy.integrate.cumtrapz to a dask array returns a numpy array, whereas the duck_array_ops.cumulative_trapezoid that I copied from duck_array_ops.trapz does return a dask array (so I assume it must not force a compute?). Maybe @fujiisoup has thoughts though?

BTW the test failures are in test_cftimeindex.py so I assume they're unrelated to these changes.

Tests pass! Think this is ready for review now @pydata/xarray

Seems like the way I've implemented surface plots requires a slightly newer version than the oldest supported matplotlib to use the "3d" projection. matplotlib-3.2.0 basically works, but some of the unit tests require matpltolib-3.3.0 to pass. I'll try to add a check and skip some tests where necessary.

Note, this PR includes the bugfix in #5099, because it needed to build on that change. #5099 is smaller, so I guess it will naturally be merged first.

I started working on this, and it seems to be almost trivial, except that plot_surface() requires it's x and y arguments to be 2d arrays. My guess is that I want to special-case the broadcasting in newplotfunc() here https://github.com/pydata/xarray/blob/ddc352faa6de91f266a1749773d08ae8d6f09683/xarray/plot/plot.py#L678-L684 but it's not immediately obvious to me how this bit works: * Why is it only ever xval that needs broadcasting? Couldn't the function have been called with a 2d x and a 1d y? * Why is it OK to check if xval.shape[0] == yval.shape[0]? I'm probably missing something, but don't see how this check can work without referring to the actual dimensions of x and y - what if the single dimension of x was actually the second dimension of y, but happened to have the same size as the first dimension of y? My first thought was that I just want to do x = x.broadcast_like(z) y = y.broadcast_like(z) actually, I think this is a bug. I'll make a new issue...

Thanks @mathause - I've merged master now.

I had the same problem. I found that da.plot.pcolormesh('lon', 'lat', shading='gouraud', infer_intervals=False) does work. Would be nice if shading='gouraud' worked by default though.

I guess the problem is that 'gouraud' needs the coordinate positions at the grid points, and I assume infer_intervals=True is creating coordinate values at the grid-cell corners (so an n x m DataArray is plotted using coordinates with lengh (n+1) and (m+1). I haven't looked into what infer_intervals actually does though. Would it be safe to just turn it off if shading='gouraud' is used?

:+1: I ran into this issue today, and found the error message strange and confusing.

I think I was trying to fix this in #4108

:+1: On a quick glance it looks to me like #4108 would fix my issue here.

I haven't managed to get set_index to line up the results. This: ``` import numpy as np import xarray as xr

ds = xr.Dataset() ds['a'] = ('x', np.linspace(0,1)) ds['b'] = ('x', np.linspace(3,4)) ds = ds.rename(b='x') ds = ds.set_coords('x') ds = ds.set_index(x='x')

print(ds) print('indexes', ds.indexes) `` Produces a Dataset that now has an indexx, but nox` coordinate. The 'assignment of 1D variable' version produces both a coordinate and an index.

One suggested solution in #2461 was to use swap_dims(), and the following does produce the desired result ``` import numpy as np import xarray as xr

ds = xr.Dataset() ds['a'] = ('x', np.linspace(0,1)) ds['b'] = ('x', np.linspace(3,4)) ds = ds.swap_dims({'x': 'b'}) ds = ds.rename(b='x')

print(ds) print('indexes', ds.indexes) `` But (1) having to useswap_dims()to create a coordinate seems bizarre, (2) I think it's a bit ugly to get rid of thexdimension withswap_dims(), and then have to rename the new dimension back tox, when what I wanted was to add a coordinate tox`.

I found the behaviour confusing because I wasn't aware of the index variables at all... I create a coordinate with set_coords(); then do a bunch of manipulations involving slicing pieces out of the Dataset and re-combining them with concat (at least I think this is the relevant part of my code...); finally test the result by taking slices again (with isel()) of the result and comparing them to slices of the original Dataset with xarray.testing.assert_identical(), which failed because of the missing indexes. I guess somewhere in the manipulations I did, some operation created the coordinate in a new Dataset object by assignment, at which point it generated an index for the coordinate too.

I may well be missing other context that makes this an undesirable thing to do, but for my use-case at least, I think it would make more sense if set_coords() created an index if the coordinate is a dimension coordinate (or whatever the actual criterion is for assignment of a 1d variable to create an index).

Apparently it's possible to work around this by using merge() instead of assigning. The following produces the expected output: ``` import numpy as np import xarray as xr

ds = xr.Dataset() ds["a"] = xr.DataArray(np.linspace(0., 1.), dims="x")

ds["x"] = xr.DataArray(np.linspace(0., 2., len(ds["x"])), dims="x") ds["x"].attrs["foo"] = "bar"

print(ds["x"])

ds["b"] = xr.DataArray(np.linspace(0., 1.), dims="x")

ds = ds.merge(xr.DataArray(np.linspace(0., 1.), dims="x", name="b"))

print(ds["x"]) ```

Unfortunately I don't have time at the moment to dig into why this is happening, but it was confusing because I found a bug where attrs had disappeared, and eventually chased it down to where some new variables were added to my Dataset, which I would not expect to change existing coordinates!

Thanks @dcherian :smile: You're very welcome!

Thanks @keewis! I think we've addressed all the review comments now.

calling on a array with missing values raises ValueError: All-NaN slice encountered (try ds = xr.tutorial.open_dataset("rasm"); ds.argmin(dim="y"))

Some missing values should be OK. In your example though for example In [27]: ds.Tair.isel(time=0, x=0).values Out[27]: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]) So there are no values in the array slice that argmin should be applied over, which is an error.

I guess we could add some special handling for this (not sure what though, because we can't set a variable with type int to nan or None), but I think that is a separate issue that would need a new PR.

ds.argmin(dim=...) returns a single index (the same result as ds.argmin()) but ds.Tair.argmin(dim=...) returns something different from ds.Tair.argmin(). Is that intentional?

I think this is a bug. dim=... to Dataset.argmin should be an error because dim=... returns a dict of results, and it's not clear how to do that consistently for a Dataset that might have several members with different combinations of dimensions.

Merge conflicts fixed, this PR should be ready to review/merge.

Yes, should be simple to correct. I have a unit test now that reproduces the error.

@pydata/xarray a question: for combine_by_coords should the default be combine_attrs="no_conflicts" (current code) or combine_attrs="drop" (current docs)?

I guess this would happen if somebody calls np.argmin(xarray_object)? Does this happen from inside xarray somewhere, or in external code?

If it's the former, then we should just fix xarray not to do this. If it's only in external code, I would still consider breaking this behavior.

test_units.py has a couple of tests that apply numpy functions to xr.Variables, xr.DataArrays and xr.Datasets. Those tests break for argmin and argmax without the out argument. Should I just remove the argmin and argmax function calls? The tests would still call them as methods.

@pydata/xarray - I think this PR is ready to be merged, are there any changes I should make? Thanks!

I eventually found that the cause of the errors I was getting was that the argmin and argmax methods did not have an out argument. In order for the methods to be wrapped by numpy (and I guess dask is the same), the call signature of np.argmin() has to be supported, which means axis and out arguments are needed.

@shoyer I've kept the use of injected functions, because removing the injected argmin and argmax completely meant re-implementing handling of skipna and only within _unravel_argminmax, which seemed less nice to me. If/when method injection is refactored, it would be nice to include a mechanism to override the core operation with an explicitly implemented, extended version like argmin/argmax.

I think this PR is ready for review now :smile:

I think indexed_array = val_arr.isel(z=z_indices) should work, if I'm understanding what you want. I haven't tested though...

Sorry for the noise, but at least I will be glad of something to search for next time I forget how to do something like this!

Thanks @lanougue, @dcherian I think I see the simple answer now: use a deepcopy first, to leave the Dataset as is, then del the DataArray when finished with it, e.g. ``` da1 = ds["variable1"].copy(deep=True) ... do stuff with da1 ... del da1

da2 = ds["variable2"].copy(deep=True) ... do stuff with da2 ... del da2

... etc. ```

Thanks @lanougue, but what if I might want to re-load da1 again later, ie. if da came from some Dataset da = ds["variable"], I want to leave ds in the same state as just after I did ds = open_dataset(...)? Wouldn't del da remove "variable" from ds? Or maybe not free any memory if ds still has a reference to the DataArray?

OK, I think I've answered my own question. Looks like dask can handle this workflow already, something like: ```

do_some_work does not call .load() or .compute() anywhere

result1a, result1b, result1c = dask.compute(do_some_work(ds["variable1"])

result2a, result2b, result2c = dask.compute(do_some_work(ds["variable2"])

... etc. ```

I do still wonder if there might be any case where .release() might be useful...

These test failures seem to have uncovered a larger issue: overriding a method injected by ops.inject_all_ops_and_reduce_methods(DataArray, priority=60) is tricky. I think there may be a way around it by mangling the name of the injected operation if that method is already defined on the class; I tried this and it nearly works, but I think I need to implement an override for argmin/argmax on Variable, because DataArray.__array_wrap__ uses the Variable version of the method... Work in progress!

The rename to _argmin_base is also causing test_aggregation (which is a test of units with pint) to fail, because the numpy function and the Variable method do not have the same name. Any suggestions how to fix this? It seems tricky because I can't pass either argmin (no Variable method) or _argmin_base (no numpy method) to the test.

I've updated so the new functionality is provided by argmin() and argmax(), when they are passed a multiple dimensions. * I've added Dataset.argmin() and Dataset.argmax(), which only work for a single dimension, but give informative exceptions for unsupported cases. If anyone has some behaviour that they would like for Dataset.argmin() with multiple dimensions, I'm happy to implement it, but returning something like a dict of Datasets didn't seem nice, and as far as I understand a Dataset cannot contain something like a dict of DataArrays. * I renamed argmin and argmax in ops.py and duck_array_ops.py to _argmin_base and _argmax_base so that I could still use them inside the argmin and argmax methods. Any alternatives, or better suggestions for naming them? * Variable no longer has argmin or argmax methods. Is that a problem? I just updated test_variable.py and test_dask.py to use _argmin_base and _argmax_base.

@shoyer I think your last option sounds good. Questions: * What should da.argmin() with no arguments do? * Currently returns the flattened index of the global minimum. * I think returning a dict of indices would be much more useful, but it does change existing behaviour (more useful because you can then do da.isel(da.argmin())). * Could anyway do da.argmin(da.dims) to get the dict result, but that's a bit ugly * Could have something like da.argmin(...) - maybe as a temporary workaround while we deprecate the current behaviour of da.argmin()? * If we overload argmin, what's the cleanest way to get the existing argmin() method within the new one? Would something like from .duck_array_ops import argmin as argmin_1d work?

Maybe worth noting, at the moment if you try to call argmin(dim=("x", "y")) with multiple dimensions, there's a not-very-helpful exception ```

array.argmin(dim=("x", "y")) Traceback (most recent call last): File "/home/john/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py", line 61, in _wrapfunc return bound(args, *kwds) TypeError: 'tuple' object cannot be interpreted as an integer

During handling of the above exception, another exception occurred:

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/<...>/xarray/xarray/core/common.py", line 46, in wrapped_func return self.reduce(func, dim, axis, skipna=skipna, kwargs) File "/<...>/xarray/xarray/core/dataarray.py", line 2288, in reduce var = self.variable.reduce(func, dim, axis, keep_attrs, keepdims, kwargs) File "/<...>/xarray/xarray/core/variable.py", line 1579, in reduce data = func(input_data, axis=axis, kwargs) File "/<...>/xarray/xarray/core/duck_array_ops.py", line 304, in f return func(values, axis=axis, kwargs) File "/<...>/xarray/xarray/core/duck_array_ops.py", line 47, in f return wrapped(args, kwargs) File "<array_function internals>", line 6, in argmin File "/<...>/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py", line 1267, in argmin return _wrapfunc(a, 'argmin', axis=axis, out=out) File "/<...>/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py", line 70, in _wrapfunc return _wrapit(obj, method, args, kwds) File "/<...>/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py", line 47, in _wrapit result = getattr(asarray(obj), method)(*args, kwds) TypeError: 'tuple' object cannot be interpreted as an integer ```

Any plans to finish/merge this PR? indexes_max and indexes_min would be very nice to have! (See also #3160.) Although I guess the idxmin/idxmax are superseded by #3871?

The Linux-py38-upstream-dev test is failing in cftime tests at various places similar to this ``` with create_tmp_file() as tmp_file: original.to_netcdf(tmp_file) with pytest.warns(None) as record: with open_dataset(tmp_file, use_cftime=False) as ds: assert_identical(expected_x, ds.x) assert_identical(expected_time, ds.time)

          assert not record

E assert not WarningsChecker(record=True) ``` Don't think it's related to this PR...

Doesn't look like it, but it seems like a thing that could be useful, e.g. for handling combining attrs dicts, so why delete it?

For specifying which object, one possibility would be to pass an int to combine_attrs in merge(), concat() or combine_by_coords(), or a tuple of int to combine_nested, giving the index of the object to use attributes from. This feature would need new tests writing though, so I'd suggest implementing it in a new PR.

Should there perhaps be another option to specify which object to get the attrs from? I'm just thinking by analogy to how open_mfdataset now lets you specify which file you want the attrs from.

That would actually be nice to have in concat() for a use-case I have. It's not immediately obvious to me how to implement it though. For merge() or concat() you could give an integer index. I think open_mfdataset used the file-name (?), but there's no equivalent for combine_by_coords or combine_nested is there?

I've added an extra test that checks the right thing happens when we do not change the bad, under or over colors. That test needs to set vmin and vmax to get sensible defaults for the under and over colors: without vmin or vmax, mpl.colors.from_levels_and_colors() is called with extend = 'neither'; then it sets _rgba_under and _rgba_over in its result to (0.0, 0.0, 0.0, 0.0) (but these values can never be used, because with no vmin or vmax there can be no data outside of the plotted range).

I think this is ready now. The unit test failures on my local computer appear to be unrelated to this PR (see earlier references to #3673 by @keewis and #3747 by @mathause).

Not sure why some tests are failing... The ones I can see are to do with CFTime. I can't see a relation to the changes in this PR??

Looks like that was it, thanks @dcherian!

Might be somethnig to do with https://github.com/pydata/xarray/blob/ed05f9862622b00f40f7b9b99ccdb0ab3766ff0f/xarray/plot/dataset_plot.py#L137 ?

The docstring for matplotlib.colors.Normalize says

If clip is True, masked values are set to 1; otherwise they remain masked. Clipping silently defeats the purpose of setting the over, under, and masked colors in the colormap, so it is likely to lead to surprises; therefore the default is clip = False.