xarray对日数据逐年份计算月均值并保留原时间维度的方法问询

我看过很多跨多年日数据计算气候态月均值的教程，但我需要对xarray中的日数据单独按每一年计算各月均值，最终得到每个经纬度格点2020年1月、2020年2月……2024年12月的对应均值。
我的xarray维度为Frozen({'time': 1827, 'lon': 180, 'lat': 90})，我尝试过代码：
var_resampled = var_diff.resample(time='1M').mean()
但它计算的是所有年份相同月份的均值，也就是2020-2024年所有1月的平均。
之后我又尝试了以下代码：

def mon_mean(x):
    return x.groupby('time.month').mean('time')

# 按年份分组，再应用上述函数
var_diff_mon = var_diff.groupby('time.year').apply(mon_mean)

这段代码的计算结果符合需求，但输出的维度变成了"month"和"year"，不再是原来的"time"维度。
请问是否有其他方法可以单独按年份计算月均值，或是有办法让上述groupby代码运行后仍保留原time维度，仅按年月聚合？
补充：我也尝试过cdo monmean指令，据我了解它同样输出的是跨多年的月平均结果。

可先使用上述groupby代码完成计算，再调用stack方法即可恢复原time维度，代码如下：

var_diff_mon.stack(time=("year", "month"))

内容的提问来源于stack exchange，提问作者random__human