How to group by actual winter season

[laurabaker-eng]

I recently discovered the function

field.collapse('T: mean',group=cf.djf())

However, I've discovered that rather than grouping by actual winter season, I think what it's doing is grouping by the Jan, Feb and Dec of the same calendar year.  Or at least, what it's producing is definitely not right!
Is there a way of doing this so that it actually gets the D,J,F from the same season?
What I've always done in the past is something along the lines of:

field_dec = field.subspace(time = cf.month(12))
field_jan = field.subspace(time = cf.month(1))
field_feb = field.subspace(time = cf.month(2))

field_djf = field_dec.copy()
field_djf.data[:-1] = (field_dec.data[:-1]+field_jan.data[1:]+field_feb.data[1:])/3.

But I'm wondering if there is a more efficient way of doing this, without manipulating the data itself?

[sadielbartholomew]

Hi @laurabaker-eng.

However, I've discovered that rather than grouping by actual winter season, I think what it's doing is grouping by the Jan, Feb and Dec of the same calendar year.

It should group by actual winter season, not the latter behaviour (the same calendar year's Dec, Jan and Feb), which would be a bug. However, I've done some testing and it seems like it does do the actual winter season collapse as we intend and document, at least for our latest version of cf-python (3.19.0), so I can't recreate the incorrect behaviour you are reporting. See sections below for my evidence.

What makes you think that the output is not right? That might help us get to the bottom of what if anything is wrong. Are weights involved in any way, since this might make the result seem different to what you expect at first?

Please can you confirm which versions of cf and its dependencies you are using by sharing the output of a cf.environment() call. If you could also share field.dump() and field_djf.dump() after the field_djf = field.collapse('T: mean',group=cf.djf()) collapse that would be really useful to see, to help us get to the bottom of this. Thanks.

Evidence that the DJF collapse is working for actual winter for v3.19.0

Method

To visually check the behaviour, I adapted one of our recipes to hack monthly time series data so the data array itself is in a 'staircase' form whereby each value is equal to the final digit of the year the data is registered for on the time axis. That way, we can see if the December that is included is that from the previous calendar year as the Jan and Feb, or the same year - because if it is from the same calendar year it will have the same value as for Jan and Feb and therefore the mean collapse will be equal to the final digit of the year. Whereas with the intended and correct behaviour, it will get the December value from the previous year and therefore be 2/3 of the way between the two digit-derived values question. (The plot and code below might make this clearer.)

Code and output

Code

import matplotlib.pyplot as plt
import cfplot as cfp

import cf

f = cf.read("~/recipes/ERA5_monthly_averaged_SST.nc")
sst = f[0]  # this gives the sea surface temperature (SST)
sst_max = sst.collapse("area: maximum")
sst_max.squeeze(inplace=True) #TODO clip

# Use this field as en example
print(sst_max)

# Hack the data so that it is a staircase like value equal to the final
# digit of the year the data is registered for on the time axis
year_final_digits = []
for dt_value in sst_max.construct("long_name=time"):
    year = [dt.year for dt in dt_value.datetime_array]
    # Extract from one-item list and take the final digit
    year_final_digit = int(str(year[0])[-1])
    year_final_digits.append(year_final_digit)

sst_max.set_data(year_final_digits)
print(sst_max.data.array)

# Do the collapse
sst_max_collapse = sst_max.collapse('T: mean', group=cf.djf())
print(sst_max_collapse.data.array)

# Plot to illustrate
cfp.gopen()
cfp.lineplot(
    sst_max,
    color="grey",
    label="All months",
)
cfp.lineplot(
    sst_max_collapse,
    color="red",
    markeredgecolor="red",
    marker="o",
    label="DJF group collapse",
    title="Testing the DJF/winter seasonal collapse",
)
cfp.gclose()

Terminal output (STDOUT)

$ python tester-lb-script-sadie.py
Field: long_name=Sea surface temperature (ncvar%sst)
----------------------------------------------------
Data            : long_name=Sea surface temperature(long_name=time(996)) K
Cell methods    : area: maximum
Dimension coords: long_name=time(996) = [1940-01-01 00:00:00, ..., 2022-12-01 00:00:00] gregorian
                : long_name=latitude(1) = [0.0] degrees_north
                : long_name=longitude(1) = [179.875] degrees_east
[0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3
 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 6 6
 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9
 9 9 9 9 9 9 9 9 9 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2
 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5
 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8
 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1
 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7
 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 0 0 0 0 0 0 0 0 0 0
 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3
 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6
 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 0
 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3
 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6
 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9
 9 9 9 9 9 9 9 9 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2
 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5
 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8
 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4
 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7
 7 7 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 0 0 0 0 0 0 0 0 0 0 0
 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3
 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 7
 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9 9 0 0
 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2]
[0.         0.66666667 1.66666667 2.66666667 3.66666667 4.66666667
 5.66666667 6.66666667 7.66666667 8.66666667 3.         0.66666667
 1.66666667 2.66666667 3.66666667 4.66666667 5.66666667 6.66666667
 7.66666667 8.66666667 3.         0.66666667 1.66666667 2.66666667
 3.66666667 4.66666667 5.66666667 6.66666667 7.66666667 8.66666667
 3.         0.66666667 1.66666667 2.66666667 3.66666667 4.66666667
 5.66666667 6.66666667 7.66666667 8.66666667 3.         0.66666667
 1.66666667 2.66666667 3.66666667 4.66666667 5.66666667 6.66666667
 7.66666667 8.66666667 3.         0.66666667 1.66666667 2.66666667
 3.66666667 4.66666667 5.66666667 6.66666667 7.66666667 8.66666667
 3.         0.66666667 1.66666667 2.66666667 3.66666667 4.66666667
 5.66666667 6.66666667 7.66666667 8.66666667 3.         0.66666667
 1.66666667 2.66666667 3.66666667 4.66666667 5.66666667 6.66666667
 7.66666667 8.66666667 3.         0.66666667 1.66666667 2.        ]

Plot output