How to calculate seasonal temporal correlations?

[krikru]

If I have two iris cubes, A and B, both on the same grid and with the same time points, how can I calculate the intra-seasonal (according the seasons that are added when calling add_season) Pearson correlations between A and B for each grid point? Can I do this somehow by only using Iris functions (i.e., not extracting the data property from both cubes and using NumPy), and without using the definition of the Pearson correlation?

Pearson correlation can usually be calculated with pearsonr, but this function doesn't seem to enable aggregation by an auxiliary variable, and I need to aggregate by season.

[rcomer]

I have a long-standing aspiration to do something like this, and a couple of failed attempts in old Iris branches. It might be easier to implement now that we have Resolve to get the data lined up properly, but I haven’t given it any deep thought recently.

As of now, the functionality doesn’t exist directly in Iris. You could write a loop to extract each season, do the correlation and add the correlation cube for each season to a CubeList. Then merge them all back together.

[krikru]

I created a workaround:

################################################################################
# IMPORTS
################################################################################


import iris
import iris.coord_categorisation
import iris.analysis.stats as ias


################################################################################
# FUNCTIONS
################################################################################


def super_MEAN(cube, coords_collapsed, coords_aggregated_by, weights=None):
    mean = cube
    if coords_collapsed:
        mean = mean.collapsed(
            coords_collapsed,
            iris.analysis.MEAN,
            weights=weights)
    if coords_aggregated_by:
        mean = mean.aggregated_by(
            coords_aggregated_by,
            iris.analysis.MEAN,
            weights=weights)
    return mean


def super_VARIANCE(cube, coords_collapsed, coords_aggregated_by, mean=None, ddof=0, weights=None):
    if mean is None:
        mean = super_MEAN(
            cube=cube,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            weights=weights)
    square_of_mean = mean**2
    mean_of_squares = super_MEAN(
        cube=cube**2,
        coords_collapsed=coords_collapsed,
        coords_aggregated_by=coords_aggregated_by,
        weights=weights)
    if ddof == 0:
        return mean_of_squares - square_of_mean
    else:
        # This requires weighted aggregation with the aggregator iris.analysis.COUNT
        raise NotImplementedError("Support for non-zero ddof in super_VARIANCE calculation not implemented yet")


def super_STD_DEV(cube, coords_collapsed, coords_aggregated_by, var=None, mean=None, ddof=0, weights=None):
    if var is None:
        var = super_VARIANCE(
            cube=cube,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            mean=mean,
            weights=weights)
    return var**0.5


def super_covariance(cube_a, cube_b, coords_collapsed, coords_aggregated_by, mean_a=None, mean_b=None, ddof=0, weights=None):
    if mean_a is None:
        mean_a = super_MEAN(
            cube=cube_a,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            weights=weights)
    if mean_b is None:
        mean_b = super_MEAN(
            cube=cube_b,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            weights=weights)
    if ddof == 0:
        mean_ab = super_MEAN(
            cube=cube_a * cube_b,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            weights=weights)
        return mean_ab - mean_a * mean_b
    else:
        # This requires weighted aggregation with the aggregator iris.analysis.COUNT
        raise NotImplementedError("Support for non-zero ddof in super_covariance calculation not implemented yet")


def super_pearsonr(cube_a, cube_b, coords_collapsed, coords_aggregated_by, mean_a=None, mean_b=None, var_a=None, var_b=None, std_dev_a=None, std_dev_b=None, cov=None, ddof=0, weights=None):
    if mean_a is None and (cov is None or std_dev_a is None):
        mean_a = super_MEAN(
            cube=cube_a,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            weights=weights)
    if mean_b is None and (cov is None or std_dev_b is None):
        mean_b = super_MEAN(
            cube=cube_b,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            weights=weights)
    if cov is None:
        cov = super_covariance(
            cube_a=cube_a,
            cube_b=cube_b,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            mean_a=mean_a,
            mean_b=mean_b,
            ddof=ddof,
            weights=weights)
    if std_dev_a is None:
        std_dev_a = super_STD_DEV(
            cube=cube_a,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            mean=mean_a,
            var=var_a,
            weights=weights)
    if std_dev_b is None:
        std_dev_b = super_STD_DEV(
            cube=cube_b,
            coords_collapsed=coords_collapsed,
            coords_aggregated_by=coords_aggregated_by,
            mean=mean_b,
            var=var_b,
            weights=weights)
    return cov / (std_dev_a * std_dev_b)

The functions take both coordinates for collapsing the data along and coordinates for aggregating the data by; hence the prefix "super_." This makes all of them, except super_MEAN, more powerful compared to only being able to calculate the aggregators by collapsing along dimensions or by aggregating by dimensions but not being able to do both simultaneously.

Is this functionality something we would like to add to Iris somehow?

[krikru]

By the way, is there some way to quickly check whether two cubes are the same (i.e., refer to the same data)? In that case, super_VARIANCE could simply be defined as a function call to super_covariance with cube provided as both cube_a and cube_b, and mean provided as both mean_a and mean_b, without loosing more than a negligible amount of time.

[acchamber]

No, there's been some thinking on if you can make cubes share data as you are suggesting, but it's not how iris is currently designed.

[trexfeathers]

#3172

[trexfeathers]

@SciTools/peloton:

Hi @krikru, thanks so much for sharing your code. At first glance, this seems like something that would be good to have in Iris API. But as you've probably guessed, that's a long list! So thank you very much for sharing your code, so that everyone can already benefit, before any Iris development work.

[bjlittle]

@SciTools/peloton This discussion has been captured in #5541

Thanks @krikru 🍻