annual mean in cmip6-ng is days-in-month-weighted

[mathause]

🤦‍♂

We assume that the mean of the monthly data equals the mean of the 12 monthly values (or that subtracting the upsampled annual data from the monthly yields anomalies around 0):

mesmer/mesmer/stats/_harmonic_model.py

Lines 352 to 353 in 4c00906

	# subtract annual mean to have seasonal anomalies around 0
	seasonal_deviations = monthly_target - yearly_predictor

However, cmip6-ng calculates the annual mean using cdo's yearmonmean which weights the data by the number of days per month. While this is more correct it violates our assumption.

Example (IPSL-CM6A-LR, 2014, full resolution):

Details

from pathlib import Path

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import mplotutils as mpl
import xarray as xr

import mesmer

# =====

model = "IPSL-CM6A-LR"
cmip6_data_path = Path("/net/atmos/data/cmip6-ng/")


path_tas_ann = cmip6_data_path / "tas" / "ann" / "g025" 
fN_hist_ann = path_tas_ann / f"tas_ann_{model}_historical_r1i1p1f1_g025.nc"

path_tas_mon = cmip6_data_path / "tas" / "mon" / "g025"
fN_hist_mon = path_tas_mon / f"tas_mon_{model}_historical_r1i1p1f1_g025.nc"

time_coder = xr.coders.CFDatetimeCoder(use_cftime=True)
tas_ann = xr.open_dataset(fN_hist_ann, decode_times=time_coder).load()
tas_mon = xr.open_dataset(fN_hist_mon, decode_times=time_coder).load()


# =====

f, ax = plt.subplots(subplot_kw={"projection": ccrs.Robinson()})


h = (tas_mon.resample(time="YE").mean().isel(time=-1).tas - tas_ann.isel(time=-1).tas).plot(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, center=0)
ax.coastlines()

mpl.colorbar(h, ax)

# ax.set_title(r"$\Delta\mathrm{T}^{\mathrm{ann}}$T ann from averaged mon minus weighted mon")
ax.set_title(r"$\Delta\mathrm{T}^{\mathrm{ann}}$ from averaged $\mathrm{T}^{\mathrm{mon}}$ $-$ weighted $\mathrm{T}^{\mathrm{mon}}$")


plt.savefig("deltaT.png", bbox_inches="tight")

[mathause]

Of course this has an influence on the coefficients (computed on the coarse grid):

but I don't know what to do about this 🙁

[mathause]

To my surprise the RSS is about 20% larger when calculating the harmonic model with the weighted averaged monthly data. (At least for the one-ensemble member coarse grid example dataset)

edit: indicated the wrong data

[mathause]

to reproduce:

# calculate annual mean from monthly
tas_y_from_m = tas_m.resample(time="YE").mean()
tas_y_from_m = tas_y_from_m.reindex_like(tas_y, method="nearest")

# stack
tas_stacked_y_from_m = mask_and_stack(tas_y_from_m, threshold_land=THRESHOLD_LAND)

# fit harmonic model
harmonic_model_fit_from_m = mesmer.stats.fit_harmonic_model(
    tas_stacked_y_from_m.tas, tas_stacked_m.tas
)

[mathause]

I also checked it with the simulations from #572 and here we find ever-so-slightly smaller rss when using the averaged annual mean

what	weighed	averaged
indiv	2503322	2500707
ave	4147570	4179860
pooled	2585725	2585043

what	weighed	averaged
historical	1542897	1540993
ssp126	335461	334997
ssp585	624964	624717