Rafael Schouten

This needs a lot of work so type promotion generally work on dimensions.

Hard to say not knowing what the errors are or what is in the vector

For RasterSeries you need to pass some dimensions that those files are over as the second argument

That sounds good, no reason it can't work. You can group by multiple dims now but as separate axes, not as a mask. I guess we should allow grouping by...

I just realised we already have that behaviour of passing a DimArray a the same as `dims(A)`. But maybe we can change that in the next breaking change to work...

That is a cool idea. How do we store the group indices for masks that big? That array could also get large. Or is it also doing online stats? (Ok...

Ok guess we will need a groupby/combine function for that. `xarray_reduce` is a pretty weird name (with no mention of the grouping part). Any idea what to call it?

I'm going to add `combine` anyway for `combine(mean, groupedarray)`, maybe we can use `combine(mean, GroupBy(A, args...; kw...))` as a lazy version? (Where `GroupBy` has identical syntax to `groupby`)

I just realised allowing this will let us have the DataFrames.jl `groupby` syntax for `AbstractDimStack` like `groupby(st, :layername)` and it will group the whole stack by that layer.

Base has custom methods for sum to do Kahan summation? We could easily add that here - simply sum each chunk separately then sum them all. Should be only minutely...