Efficient cube multiple statistic calculations

[bjlittle]

Extend the capability of iris to allow the user to efficiently calculate multiple statistics from the same source cube when using iris.cube.Cube.collapsed and iris.cube.Cube.aggregated_by.

For example:

from iris.analysis import MEAN, SUM, MAX
cubes = cube.collapsed('model_level_number', [MEAN, SUM, MAX])

As opposed to the less efficient:

cubes = [cube.collapsed('model_level_number', agg) for agg in [MEAN, SUM, MAX]]

Ditto for iris.cube.Cube.aggregated_by.

The resolution of source cube data is only on the increase. As such, iris needs to be sympathetic to those scientists that require to calculate multiple statistics and analyses over high and very high resolution data by offering an efficient workflow of streaming once, calculating many. As opposed to streaming once, calculating once.

[bjlittle]

@BenMGeo do you have some concrete user examples from ESMValTool that we could reference here?

[bjlittle]

This might also be an opportunity for us to rationalise and unify the code base of collapsed and aggregated_by.... we've been carrying that technical debt for way, way too long 😉

[BenMGeo]

The simplest example coming to my mind is calculating trends along time axis (or pressure levels).
Usually, the implementations of trends have more results than just the slope information (first priority), but also the p-values, R², correlation, ...
It would be gold to retrieve all these cubes at once. Currently, I wrote a wrapper where I request one single result after the other, which takes a lot of time.

I can find and drop code snippets from my implementation if requested.

[pp-mo]

I think this is now at least partially addressed by #2967 / #3013
I.E. the CubeList.realise_data function

So, instead of the above inefficient example

cubes = [cube.collapsed('model_level_number', agg) for agg in [MEAN, SUM, MAX]]

where the inefficient bit is really when you ...

for cube in cubes:
    cube.data

You can now write ...

cubes = [cube.collapsed('model_level_number', agg) for agg in [MEAN, SUM, MAX]]
CubeList(cubes).realise_data()

[pp-mo]

RE: my previous.
We made this into a CubeList method, which was not my personal preference, but in any case I was always a bit doubtful about how to publicise this.
Any better ideas ?

[rcomer]

I may have misunderstood, but it looks to me like @BenMGeo's example is a bit different from the OP. I think the trend analysis is a single numpy/scipy function which returns several arrays. So we'd want a single aggregator to wrap that function and return a list of cubes?

[pp-mo]

@rcomer we'd want a single aggregator to wrap that function and return a list of cubes?

In that case it will require a dask-specific solution

@BenMGeo I can find and drop code snippets from my implementation if requested.

👍 👍 😁

[BenMGeo]

@rcomer yes, that's right.
I'd imagine something like that:

def something(array, **kwargs):
  a = suba(array, **kwargs)
  b = subb(array, **kwargs)
  return a, b

C = iris.cube.Cube(...)

A, B = C.collapsed("time", wrapper(something()), **kwargs) 

Not sure if you understand this "pseudocode".

[BenMGeo]

I'm using it for trends here, calculating slope and p-value:

{cube is iris.cube.Cube and has dimensions ["time","latitude","longitude"]}

...
        # simple linear trend (slope) and p-values
_, S, _, P = utils.temporal_trend(cube, pthres=1.01)
...

referring to a submodule utils:

def temporal_trend(cube, pthres=1.01):
    # time difference in days and decades (not generel solution)
    dt = np.asarray(cube.coord('time').points - cube.coord('time').points[0]
                    ).astype('float') / 365.2425 / 10
    # ensure that a masked array is used
    x = np.ma.array(dt, mask=dt != dt)

    R, S, I, P, C = __corr_single__(cube, x, pthres=pthres)

    R.long_name = cube.long_name + ' (correlation)'
    S.long_name = cube.long_name + ' ($\partial x / \partial t$)'
    I.long_name = cube.long_name + ' (offset)'
    P.long_name = cube.long_name + ' (p-value)'

    if S.units in [None, 'no_unit', '1', 'unknown']:
        S.units = cf_units.Unit('0.1 year-1')
    else:
        S.units = cf_units.Unit(str(S.units) + ' 0.1 year-1')
    R.units = '-'
    I.units = cube.units
    P.units = 1
    return R, S, I, P

Currently, in __corr_single__ I calculate the data based on the core_data with the method stats.mstats.linregress along a spatial loop which returns r, s, i and p.

It might be an issue to align the units by default, but unit-less return of a list of cubes with further notice/warning might be feasible (?)