I am interested in the coarsen with weights scenario that @dcherian and @mathause described here for a current project of ours.

I solved the issue manually and its not that hard ```python import xarray as xr import numpy as np

example data with weights

data = np.arange(16).reshape(4,4).astype(float)

add some nans

data[2,2] = np.nan data[1,1] = np.nan

create some simple weights

weights = np.repeat(np.array([[1,2,1,3]]).T, 4, axis=1) weights

da = xr.DataArray(data, dims=['x', 'y'], coords={'w':(['x','y'], weights)}) da ```

python masked_weights = da.w.where(~np.isnan(da)) # .weighted() already knows how to do this da_weighted = da * masked_weights da_coarse = da_weighted.coarsen(x=2, y=2).sum() / masked_weights.coarsen(x=2, y=2).sum() da_coarse

but I feel all of this is duplicating existing functionality (e.g. the masking of weights based on nans in the data) and might be sensibly streamlined into something like: python da.weighted(da.w).coarsen(...).mean() at least from a user perspective (there might be unique challenges with the implementation that I am overlooking here).

Happy to help but would definitely need some guidance on this one.

I do believe that this would provide a very useful functionality for many folks who work with curvilinear grids and want to prototype things that depend on some sort of scale reduction (coarsening).

Also cc'ing @TomNicholas who is involved in the same project 🤗

This should be fully reproducible in the pangeo cloud deployment, but is unfortunately only available as 'requester pays' for other local machines.=

How many variables are in ds. You're diagnosing the graph construction time in that %timeit statement. This will scale with number of variables.

Ahh good catch. @shanicetbailey can you try to drop all but two variables from the cloud dataset (ds[['SST', 'U']]) and check again?

yes all of the grid methods (grid.diff etc) are now internally using grid_ufuncs. The axis methods are still going through the old code path, but will be deprecated soon! Please let us know how you get along with the new functionality, we are very curious for user feedback!

Hi @miniufo et al., just my two cents:

This is simpler and do not make heavy dependence of the third-party package like xgcm.

That is a fair point, but I think there is a counterpoint to be made, that xgcm gives you some more functionality (especially with the new grid_ufuncs feature) with regard to array padding. As you note, this is not needed for your particular setup, but if you use xgcm, you would get the same functionality + at a later point you might get padding on complex grid topologies for free down the line. So in the end this seems like a tradeoff between adding more dependencies vs flexibility and generalizability in the future.

I'll give a try with differentiate() and pad() to implement grad/div/vor... But some designs in xgcm also inspire me to make things much natural.

This makes me think that you really want xgcm, because these properties will naturally be located on staggered grid positions, even if your data is originally on a A grid. And once you start to try to handle these cases it would appear to me that you duplicate some of the functionality of xgcm?

I am still worried about the metrics concept introduced by xgcm. I think this should be discussed over xgcm's repo.

I second others here and think it would be great to elaborate on this on the xgcm issue tracker. But I also want to point out, that using the metrics functionality is entirely optional in xgcm, so if you desire, you can roll your own logic on top of grid.diff/interp etc.

I am very interested in this sort of functionality as an xarray accessor. If I can help in any way, please let me know. Ideally this work would come in very handy to visualize Oxygen Minimum Zones in the global ocean as isosurfaces of a 3D oxygen array.

Strictly speaking, the values are different I guess. However I think this error would be clearer if it said that the dimension order was different but the values are equal once the dimensions are transposed.

I guess this comes down a bit to a philosophical question related to @benbovy s comment above. You can either make this operation be similar to the numpy equivalent (with some more xarray specific checks) or it can check whether the values at a certain combo of labels are the same/close.

The latter would be the way I think about data in xarray as a user. To me the removal of axis logic (via labels) is one of the biggest draws for myself, but importantly I also pitch this as one of the big reasons to switch to xarray for beginners.

I would argue that a 'strict' (numpy style) comparision is less practical in a scientific workflow and we do have the numpy implementation to achieve that functionality. So I would ultimately argue that xarray should check closeness between values at certain label positions by default.

However, this might be very opinionated on my end, and a better error message would already be a massive improvement.

That might be a bit too special to be included directly in xarray. Note, however, that instead of a name you can also pass a function to combine_attrs so you should be able to customize this.

👀

Ideally, this: ```python ds1 = xr.Dataset(attrs={'a':[5]}) ds2 = xr.Dataset(attrs={'a':5})

xr.merge([ds1, ds2], combine_attrs='drop_conflicts') ```

would actually not be dropped but resolved to either 5 or [5]?

Should I close this one then?

I just stumbled over this in the cmip6_preprocessing CI. I would really appreciate a bugfix release. Cheers.

the question is whether the chunk() method should delete existing chunks attributes from encoding.

IMO this is the user-friendly thing to do.

Just ran into this issue myself and just wanted to add a +1 to stripping the encoding when .chunk() is used.

I encountered a similar problem, which I could solve by dropping all scalar non-dim coords in my dataset (works in this particular workflow of mine, but is generally not ideal). Is the solution proposed by @chrisroat general enough to be implemented? Or is there another way to avoid this situation (besides dropping the coordintates?).

Thanks for all the help! I made the entry to whats-new.rst.

do you wan't to add a what's new entry? In any case we can merge this in a few days unless someone else has a comment.

Not sure. Would you think that this is significant enough? If yes, Id be happy to do it.

Oops. Thanks for catching that. Should be fixed now.

As @dcherian pointed out above copy(..., deep=False) does fix this for all cases I am testing.

So I have added a test in #4668 and it confirms that this behavior is only occurring if the resample interval is smaller or equal than the chunks. If the resample interval is larger than the chunks it stays completely lazy...not sure if this is a general limitation? Does anyone have more insight into how resample handles this kind of workflow?

Oh I remember that too, and I didn't understand it at all...

That was it. Sorry for the false alarm! Closing this.

Do you have a suggestion how to test this? Should I write a test involving resample + weighted?

Sweet. Ill try to apply this fix for my workflow now. Happy to submit a PR with the suggested changes to weighted.py too.

Oh nooo. So would you suggest that in addition to #4559, we should have a kwarg to completely skip this?

Seems like all the other test are passing (minus the two upstream problems discussed before).

I am not understanding why that MinimumVersionPolicy test is failing (or what it does haha). Is this something upstream, or should I fix?

Ok I think this should be good to go. I have implemented all the requested changes. The remaining failures are related to other problems upstream (I think). Anything else I should add here?

Similarly on MacOSX py38 this fails: TestDask.test_save_mfdataset_compute_false_roundtrip

I am getting some failures for py38-upstream (TestDataset.test_polyfit_warnings) that seem unrelated?

Do you think this works or are further changes needed? Many thanks for the guidance so far!

I did have to fiddle with this a bit. I did change .isnull() to np.isnan() and replace the data of weights, otherwise the error would not be raised for the dask array at all. This does not look terribly elegant to me, but passes tests locally. Waiting for the CI to see if the others pass aswell. Happy to make further changes, and thanks for all the input already.

The Ci environments without dask are failing. Should I add some pytest skip logic, or what is the best way to handle this?

Another option would be to put the check in a .map_blocks call for dask arrays. This would only run and raise at compute time.

Uh that sounds great actually. Same functionality, no triggered computation, and no intervention needed from the user. Should I try to implement this?

Thanks @mathause , I was wondering how much of a performance trade off .fillna(0) is on a dask array with no nans, compared to the check.

I favor this, since it allows slicing before the calculation is triggered: I have a current situation where I do a bunch of operations on a large multi-model dataset. The weights are time and member dependent and I am trying to save each member separately. Having the calculation triggered for the full dataset is problematic and fillna(0) avoids that (working with a hacked branch where I simply removed the check for nans)

Sounds good. I'll see if I can make some time to test and put up a PR this week.

What would happen in this case if a dask array with nans is passed? Would this somehow silently influence the results or would it not matter (in that case I wonder what the check was for). If this could lead to undetected errors I would still consider a kwargs a safer alternative, especially for new users?

Sorry if my initial issue was unclear. So you favor not having a 'skip' kwarg to just internally skipping the call to .any() if weights is a dask array?

Just tried this with the newest dask version and can confirm that I do not get huge chunks anymore IF i specify dask.config.set({"array.slicing.split_large_chunks": True}). I also needed to modify the example to exceed the internal chunk size limitation: ```python import numpy as np import xarray as xr import dask dask.config.set({"array.slicing.split_large_chunks": True})

short_time = xr.cftime_range('2000', periods=12) long_time = xr.cftime_range('2000', periods=120)

data_short = np.random.rand(len(short_time)) data_long = np.random.rand(len(long_time)) n=1000 a = xr.DataArray(data_short, dims=['time'], coords={'time':short_time}).expand_dims(a=n, b=n).chunk({'time':3}) b = xr.DataArray(data_long, dims=['time'], coords={'time':long_time}).expand_dims(a=n, b=n).chunk({'time':3})

a,b = xr.align(a,b, join = 'outer') `` with the option turned on I get this fora`;

with the defaults, I still get one giant chunk.

Ill try this soon in a real world scenario described above. Just wanted to report back here.

Wondering if this is still an issue. I dont have the data to check it but in my experience these kind of operations have been much better in recent versions. Ill close this for now.

Closing in favor of SciTools/nc-time-axis#44

Apologies for not noticing this earlier. I (wrongly) assumed that xarray would handle the axis limits. Should I raise an issue over there? I am currently quite busy so it might take a while for me to be able to work on a PR.

It looks like you can set xlim in da.plot.line - could you test this?

I tried this: ```

This needs a good amount of dask workers!

ds = data_dict['CMIP.CSIRO.ACCESS-ESM1-5.piControl.Omon.gn'] ts.plot(xlim=['0005', '2000']) ``` and am getting the same error.

I can confirm that this issue is resolved for my project. Seems to not make a difference in speed anymore whether I assign the dataarray as coordinate or data variable. Thanks for the fix!

I think this issue was actually a dupe. I remember you pointing me to changes in 14.x, that improved the performance, but I cant find the other issue right now. I will have an opportunity to test this in the coming days on some huge GFDL data

I can give it a shot if you could point me to the appropriate place, since I have never messed with the dask internals of xarray.

I believe this was fixed in a recent version. Closing

This is the chunk setup

Might this be a problem resulting from numpy.vectorize?

I am having a similar problem. This impacts some of my frequently used code to compute correlations.

Here is a simplified example that used to work with older dependencies: ``` import xarray as xr import numpy as np from scipy.stats import linregress

def _ufunc(aa,bb): out = linregress(aa,bb) return np.array([out.slope, out.intercept])

def wrapper(a, b, dim='time'): return xr.apply_ufunc( _ufunc,a,b, input_core_dims=[[dim], [dim]], output_core_dims=[["parameter"]], vectorize=True, dask="parallelized", output_dtypes=[a.dtype], output_sizes={"parameter": 2},) ```

This works when passing numpy arrays:

a = xr.DataArray(np.random.rand(3, 13, 5), dims=['x', 'time', 'y']) b = xr.DataArray(np.random.rand(3, 5, 13), dims=['x','y', 'time']) wrapper(a,b)

But when I convert both arrays to dask arrays, I get the same error as @smartass101.

wrapper(a.chunk({'x':2, 'time':-1}),b.chunk({'x':2, 'time':-1}))

This used to work like a charm...I however was sloppy in testing this functionality (a good reminder always to write tests immediately 🙄 ), and I was not able to determine a combination of dependencies that would work. I am still experimenting and will report back

Could this behaviour be a bug introduced in dask at some point (as indicated by @smartass101 above)? cc'ing @dcherian @shoyer @mrocklin

EDIT: I can confirm that it seems to be a dask issue. If I restrict my dask version to <2.0, my tests (very similar to the above example) work.

Just stumbled across this issue. Is there a recommended workaround?

I am usually doing this (specific to seasons): import xarray as xr ds = xr.tutorial.open_dataset('air_temperature') airtemp_seasonal = ds.groupby('time.season').mean('time').sortby(xr.DataArray(['DJF','MAM','JJA', 'SON'],dims=['season'])) Thought this might help some folks who need to solve this problem.

You guys are the best! Thanks.

Julius J.M. Busecke, Ph.D. (he/him)

———————————————

Postdoctoral Research Associate

Princeton University • Geosciences

408A Guyot Hall • Princeton NJ

juliusbusecke.com

On Oct 29, 2019, at 10:47 AM, Deepak Cherian notifications@github.com wrote:



Totally fixed by #3453https://github.com/pydata/xarray/pull/3453!!! Both statements take the same time on that branch.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/pydata/xarray/issues/3454?email_source=notifications&email_token=ADNGY72OB7JSHN4OALG7RALQRBEGXA5CNFSM4JGIQ3AKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECQYRBQ#issuecomment-547457158, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADNGY7ZRON5TWT4GMB4HEHLQRBEGXANCNFSM4JGIQ3AA.

I think this was closed via #3338. Closing

Glad that this orphaned test (we ended up removing it, because the function was not called anymore) was still useful!

And many thanks to @dcherian for suggesting to test xgxm with the upstream master!

_{Sent with GitHawk}

Thanks for this quick implementation @dcherian. I will work on implementing testing with the xarray master downstream so we can catch these earlier.

Wooooow. Thanks. Ill have to give this a whirl soon.

If I understand correctly then it gets piped through cmap_kwargs, which seems odd to me. Do you agree that this is a bug? Or am I missing a case where it would be preferable to pass extend directly to plot.contourf?

Also FYI I have a PR open that will enable xmovie to write movie files (by invoking ffmpeg 'under the hood'). Just wanted to mention it since this might come in handy as another export option for this feature later on.

This looks amazing! Which problem are you referring to specifically @rabernat?

Cool. Ill try to give that a try some time soon

How about this one?

I am still a bit puzzled over the ds[['x', 'temperature']] suggestion. Maybe I am not getting something, but x is a dimension and the output of ds[['x', 'temperature']] is the same as ds[['temperature']]. I think it would be clearer to add a third variable (maybe in the code block above) and then select two out of three data variables?

Oh I see. It returns None if any of the keys is not found. That might indeed lead to confusion. So should I just add an example with multiple variables using ds[['var1', 'var2']]?

Ok I have submitted a PR for the xarray.Dataset.get function.

@dcherian, I was not able to find that issue you mentioned. I would certainly be interested to have a look in the future.

Point taken. I am still not thinking general enough :-)

Are we going to require that the argument to weighted is a DataArray that shares at least one dimension with da?

This sounds good to me.

With regard to the implementation, I thought of orienting myself along the lines of groupby, rolling or resample. Or are there any concerns for this specific method?

I have to say that I am still pretty bad at thinking fully object orientented, but is this what we want in general? A subclass of xr.DataArray which gets initialized with a weight array and with some logic for nans then 'knows' about the weight count? Where would I find a good analogue for this sort of organization? In the rolling class?

I like the syntax proposed by @jhamman above, but I am wondering what happens in a slightly modified example: ```

da.shape (72, 10, 15) da.dims ('time', 'x', 'y') weights = some_func_of_x(x) da.weighted(weights).mean(dim=('x', 'y')) `` I think we should maybe build in a warning that when theweights` array does not contain both of the average dimensions?

It was mentioned that the functions on ...weighted(), would have to be mostly rewritten since the logic for a weigthed average and std differs. What other functions should be included (if any)?

Found this issue due to @rabernats blogpost. This is a much requested feature in our working group, and it would be great to build onto it in xgcm aswell. I would be very keen to help this advance.

Wow. Thats awesome. Had no clue about it. I will put in a PR for the docs for sure. Might take a bit though. Ill also take a look at the keep_vars option.

Could you think of a way I would be able to diagnose this further? Sorry for these wide questions but I am not very familiar with these xarray internals.

Oh yeah, that seems totally fair to me. Thanks for clarifying. Cant wait to have this functionality! Thanks @spencerkclark

Oh sweet, I just encountered this problem. Would this work on a large dask array with a non-dask time dimension?

I have run into this problem multiple times. The latest example I found were some [CORE ocean model runs] (https://rda.ucar.edu/datasets/ds262.0/index.html#!description). The time dimension of some (they mix units) of these files is given as netcdf MRI-A_sigma0_monthly { dimensions: level = 51 ; latitude = 368 ; longitude = 364 ; time = UNLIMITED ; // (720 currently) time_bnds = 2 ; variables: double latitude(latitude) ; latitude:units = "degrees_north " ; latitude:axis = "Y" ; double longitude(longitude) ; longitude:units = "degrees_east " ; longitude:axis = "X" ; double level(level) ; level:units = "m " ; level:axis = "Z" ; double time(time) ; time:units = "years since 1948-1-1 00:00:00 " ; time:axis = "T" ; time:bounds = "time_bnds" ; time:calendar = "noleap" ; double time_bnds(time, time_bnds) ; float sigma0(time, level, latitude, longitude) ; sigma0:units = "kg/m^3 " ; sigma0:long_name = "Monthly-mean potential density (sigma-0) " ; sigma0:missing_value = -9.99e+33f ; }

I understand that 'fully' supporting to decode this unit is hard and should probably addressed upstream.

But I think it might be useful to have a utility function that converts a dataset with these units into someting quickly useable with xarray? E.g. one could load the dataset with ds = xr.open_dataset(..., decode_times=False) and then maybe call xr.decode_funky_units(ds, units='calendaryears', ...), which could default to the first day of a year (or the first day of a month for units of months since.

This way the user is aware that something is not decoded exactly, but can work with the data. Is this something that people could see useful here? Id be happy to give an implementation a shot if there is interest.

Is there anything else that I need to do at this point? Sorry for the xarray noob question...

Awesome. Just added the line. Let me know if you think it is appropriate.

Thanks. I updated the PR accordingly.

Seems like the travis builds all pass, wohoo. Please let me know if anything else is needed.

I think I have addressed all the above remarks (Many thanks for the thorough review and tips). Waiting for the CI again.

I am smh getting these errors in the backend part of the tests:

``` ======================================================== FAILURES ========================================================= _____ TestNetCDF3ViaNetCDF4Data.test_encoding_same_dtype ______

self = <xarray.tests.test_backends.TestNetCDF3ViaNetCDF4Data object at 0xd238d9a20>

def test_encoding_same_dtype(self):
    ds = Dataset({'x': ('y', np.arange(10.0, dtype='f4'))})
    kwargs = dict(encoding={'x': {'dtype': 'f4'}})
    with self.roundtrip(ds, save_kwargs=kwargs) as actual:

      assert actual.x.encoding['dtype'] == 'f4'

E AssertionError: assert dtype('>f4') == 'f4'

xarray/tests/test_backends.py:853: AssertionError _______ TestGenericNetCDFData.test_encoding_same_dtype ________

self = <xarray.tests.test_backends.TestGenericNetCDFData object at 0xd238e0588>

def test_encoding_same_dtype(self):
    ds = Dataset({'x': ('y', np.arange(10.0, dtype='f4'))})
    kwargs = dict(encoding={'x': {'dtype': 'f4'}})
    with self.roundtrip(ds, save_kwargs=kwargs) as actual:

      assert actual.x.encoding['dtype'] == 'f4'

E AssertionError: assert dtype('>f4') == 'f4'

xarray/tests/test_backends.py:853: AssertionError ``` They do not always show up...not sure what to make of it, but could be an issue with my local environment. Lets see if the CI shows this aswell.

Ok I think I have most of the things covered. All test pass for me locally. What should I add to the whats-new.rst. I thought of something like this under Enhancements (or would this be considered a bug fix?): Internal plotting now supports cftime.datetime objects as time axis (@spencerkclark, @jbusecke)

Great idea to simplify @spencerkclark. Thanks. Regarding the tests. I have removed the following:

``` @requires_cftime def test_plot_cftime_coordinate_error(): cftime = _import_cftime() time = cftime.num2date(np.arange(5), units='days since 0001-01-01', calendar='noleap') data = DataArray(np.arange(5), coords=[time], dims=['time']) with raises_regex(TypeError, 'requires coordinates to be numeric or dates'): data.plot()

@requires_cftime def test_plot_cftime_data_error(): cftime = _import_cftime() data = cftime.num2date(np.arange(5), units='days since 0001-01-01', calendar='noleap') data = DataArray(data, coords=[np.arange(5)], dims=['x']) with raises_regex(NotImplementedError, 'cftime.datetime'): data.plot() ``` And the test suite passes locally.

But I assume Ill have to add another test dataset with a cftime.datetime time-axis, which then gets dragged through all the plotting tests? Where would I have to put that in?

Many thanks for all the help

I have quickly looked into the testing and found an oddity that might be important if nc-time-axis is not installed.

So in the definition of plot in plot.py I have changed

if contains_cftime_datetimes(darray): to if any([contains_cftime_datetimes(darray[dim]) for dim in darray.dims]): Because if I understand correctly, the previous statement only checks the dtype of the actual data, not the dimensions. Is this appropriate or am I misunderstanding the syntax? In my example above it doesnt matter, because this only spits out an error message when nc-time-axis is not available.

Ok so the plotting works now with both timeseries and 2d data as follows

``` import xarray as xr import numpy as np %matplotlib inline

Create a simple line dataarray with cftime

time = xr.cftime_range(start='2000', periods=4, freq='1H', calendar='noleap') data = np.random.rand(len(time)) da = xr.DataArray(data, coords=[('time', time)]) da.plot() ![download](https://user-images.githubusercontent.com/14314623/51777752-6d427580-20cc-11e9-9c3d-a91d31b10312.png)

Check with 2d data

time = xr.cftime_range(start='2000', periods=6, freq='2MS', calendar='noleap') data2 = np.random.rand(len(time), 4) da2 = xr.DataArray(data2, coords=[('time', time), ('other', range(4))]) da2.plot() ```

Cool. Ill give it a shot right now.

Sounds good to me.

Best

Julius On Jan 23, 2019, 12:56 PM -0500, Spencer Clark notifications@github.com, wrote:

I agree @dcherian; I just pinged the PR again, but if there is no activity there by this time next week, I think we should probably move forward here. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

Is there still interest in this PR? Or did the upstream changes move ahead? I am finding myself explaining workarounds for this too students in the department, so maybe my time would be better invested getting this fix to the full community?

But obviously if things are going to be fixed upstream soon, I would devote time to other projects. Thoughts?

Oh shoot, I now remember seeing this. If this will be implemented soon I guess the PR can be discarded. Any chance you would have a quick solution for the pcolormesh plot error (second example in the PR) @spencerkclark?

One of the more general questions I had was if we should expose the conversion using nc-time-axis in the public API. That way users could easily plot the data in matplotlib, e.g.: da_new = da.convert_cftime() plt.plot(da_new.time, da_new) Just an idea...

I have been along the lines of a short example. This works for timeseries data.

``` import xarray as xr import numpy as np %matplotlib inline

Create a simple line dataarray with cftime

time = xr.cftime_range(start='2000', periods=6, freq='2MS', calendar='noleap') data = np.random.rand(len(time)) da = xr.DataArray(data, coords=[('time', time)]) da.plot() ```

For pcolormesh plots this still fails.

```

Create a simple line dataarray with cftime

time = xr.cftime_range(start='2000', periods=6, freq='2MS', calendar='noleap') data2 = np.random.rand(len(time), 4) da2 = xr.DataArray(data2, coords=[('time', time), ('other', range(4))]) da2.plot() ```

Perhaps @spencerkclark has an idea how to deal with differencing cftime.datetime objects?

I have taken a swing at restoring the internal plotting capabilities in #2665. Feedback would be very much appreciated since I am still very unfamiliar with the xarray plotting internals.

I think mean would be a good default (thinking about cell center dimensions like longitude and latitude) but I would very much like it if other functions could be specified e. g. for grid face dimensions (where min and max would be more appropriate) and other coordinates like surface area (where sum would be the most appropriate function).

On Nov 18, 2018, at 11:13 PM, Ryan Abernathey notifications@github.com wrote:

What would the coordinates look like?

apply func also for coordinate always apply mean to coordinate If I think about my applications, I would probably always want to apply mean to dimension coordinates, but would like to be able to choose for non-dimension coordinates.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

My favorite would be da.coarsen({'lat': 2, 'lon': 2}).mean(), but all the others sound reasonable to me. Also +1 for consistency with resample/rolling/groupby.

I agree with @rabernat, and favor the index based approach. For regular lon-lat grids its easy enough to implement a weighted mean, and for irregular spaced grids and other exotic grids the coordinate based approach might lead to errors. To me the resample API above might suggest to some users that some proper regridding (a la xESMF) onto a regular lat/lon grid is performed.

‚block_reduce‘ sounds good to me and sounds appropriate for non-dask arrays. Does anyone have experience how ‚dask.coarsen‘ compares performance wise?

I should add that I would be happy to work on an implementation, but probably need a good amount of pointers.

Here is the implementation that I have been using (only works with dask.arrays at this point).

Should have posted that earlier to avoid @rabernat s zingers over here. ```python def aggregate(da, blocks, func=np.nanmean, debug=False): """ Performs efficient block averaging in one or multiple dimensions. Only works on regular grid dimensions. Parameters ---------- da : xarray DataArray (must be a dask array!) blocks : list List of tuples containing the dimension and interval to aggregate over func : function Aggregation function.Defaults to numpy.nanmean Returns ------- da_agg : xarray Data Aggregated array Examples -------- >>> from xarrayutils import aggregate >>> import numpy as np >>> import xarray as xr >>> import matplotlib.pyplot as plt >>> %matplotlib inline >>> import dask.array as da >>> x = np.arange(-10,10) >>> y = np.arange(-10,10) >>> xx,yy = np.meshgrid(x,y) >>> z = xx2-yy2 >>> a = xr.DataArray(da.from_array(z, chunks=(20, 20)), coords={'x':x,'y':y}, dims=['y','x']) >>> print a <xarray.DataArray 'array-7e422c91624f207a5f7ebac426c01769' (y: 20, x: 20)> dask.array<array-7..., shape=(20, 20), dtype=int64, chunksize=(20, 20)> Coordinates: * y (y) int64 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 * x (x) int64 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 >>> blocks = [('x',2),('y',5)] >>> a_coarse = aggregate(a,blocks,func=np.mean) >>> print a_coarse <xarray.DataArray 'array-7e422c91624f207a5f7ebac426c01769' (y: 2, x: 10)> dask.array<coarsen..., shape=(2, 10), dtype=float64, chunksize=(2, 10)> Coordinates: * y (y) int64 -10 0 * x (x) int64 -10 -8 -6 -4 -2 0 2 4 6 8 Attributes: Coarsened with: <function mean at 0x111754230> Coarsenblocks: [('x', 2), ('y', 10)] """ # Check if the input is a dask array (I might want to convert this # automaticlaly in the future) if not isinstance(da.data, Array): raise RuntimeError('data array data must be a dask array') # Check data type of blocks # TODO write test if (not all(isinstance(n[0], str) for n in blocks) or not all(isinstance(n[1], int) for n in blocks)):

    print('blocks input', str(blocks))
    raise RuntimeError("block dimension must be dtype(str), \
    e.g. ('lon',4)")

# Check if the given array has the dimension specified in blocks
try:
    block_dict = dict((da.get_axis_num(x), y) for x, y in blocks)
except ValueError:
    raise RuntimeError("'blocks' contains non matching dimension")

# Check the size of the excess in each aggregated axis
blocks = [(a[0], a[1], da.shape[da.get_axis_num(a[0])] % a[1])
          for a in blocks]

# for now default to trimming the excess
da_coarse = coarsen(func, da.data, block_dict, trim_excess=True)

# for now default to only the dims
new_coords = dict([])
# for cc in da.coords.keys():
warnings.warn("WARNING: only dimensions are carried over as coordinates")
for cc in list(da.dims):
    new_coords[cc] = da.coords[cc]
    for dd in blocks:
        if dd[0] in list(da.coords[cc].dims):
            new_coords[cc] = \
                new_coords[cc].isel(
                    **{dd[0]: slice(0, -(1 + dd[2]), dd[1])})

attrs = {'Coarsened with': str(func), 'Coarsenblocks': str(blocks)}
da_coarse = xr.DataArray(da_coarse, dims=da.dims, coords=new_coords,
                         name=da.name, attrs=attrs)
return da_coarse

```

Is this feature still being considered? A big +1 from me.

I wrote my own function to achieve this (using dask.array.coarsen), but I was planning to implement a similar functionality in xgcm, and it would be ideal if we could use an upstream implementation from xarray.

I would prefer axis_aspect to allow other than square aspects.

I like the idea. I would prefer aspect_square=True to be more clear.

I am unsure if that failure is due to a time out or changes in the PR. Is ther anything else that I should change before merge? Again, many thanks for the help in getting this forward.

Yikes. Also sorry for this merge fail commit. I am sitting in a conference :-).

Awesome. Thanks a lot for the feedback @dcherian.

Yes that is probably more consistent @dcherian. Ill change it in a bit.

It seems like the problems I encountered during testing are caused by my local setup after all. The Travis CI passed. Just added the test for the divergent colormap.

I took a shot at it in #2397. Setting the options works locally on my laptop, but I am not clear how to properly test it.

I encountered this problem right now with the xarray built-in plotting. Does anybody know a workaround for the xarray plotting by any chance?

This looks good to me! thanks for the implementation! This will save lots of seconds that add up ;-)

Awesome, thanks for the clarification. I just looked at #1981 and it seems indeed very elegant (in fact I just now used this approach to parallelize printing of movie frames!) Thanks for that!