这个解决方案使用terra而不是raster，因为根据这两个包的作者的说法，“它可以做更多的事情，更快，更容易使用”。我也会做一个最小的例子，使用lubridate来处理日期。

library('terra')
library('lubridate')

r = rast(ncols = 10, nrows = 10)
values(r) = 1:ncell(r)
s = c(r, r, r, r, r, r, r, r)
s = s * 1:8
names(s) = c(
  '2020-01-01',
  '2020-01-02',
  '2020-02-01',
  '2020-02-02',
  '2021-01-01',
  '2021-01-02',
  '2021-02-01',
  '2021-02-02'
)
s

class       : SpatRaster 
    dimensions  : 10, 10, 8  (nrow, ncol, nlyr)
    resolution  : 36, 18  (x, y)
    extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    coord. ref. : lon/lat WGS 84 
    source      : memory 
    names       : 2020-01-01, 2020-01-02, 2020-02-01, 2020-02-02, 2021-01-01, 2021-01-02, ... 
    min values  :          1,          2,          3,          4,          5,          6, ... 
    max values  :        100,        200,        300,        400,        500,        600, ...

因此s在两年中的每一年中的每一个月中的每一天中都有一个层。现在我们可以使用terra::tapp()来计算聚合堆栈，一次一个。本示例将写入内存，但对于真实的情况，最好考虑使用filename参数将这些结果写入磁盘。
一个二个一个一个
现在我们每个月都有一个层，包含每天的总和，并避免了日历问题（闰年），这要归功于lubridate。

# yearly stack (mean of every month)
# this time we have to deal with the names generated by `terra`
dts = names(monthly_s) |> gsub('X', '', x = _) |> lubridate::ym()
yearly_i = lubridate::year(dts) |> factor()

# aggregate the stack
yearly_s = terra::tapp(monthly_s, index = yearly_i, fun = mean)
yearly_s

class       : SpatRaster 
    dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
    resolution  : 36, 18  (x, y)
    extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    coord. ref. : lon/lat WGS 84 
    source      : memory 
    names       : X2020, X2021 
    min values  :     5,    13 
    max values  :   500,  1300

注意，这里的初始堆栈有漂亮的名字，这段代码应该根据需要进行调整以处理生成的堆栈。

示例数据

library(terra)
r <- rast(ncols=10, nrows=10, nlyr=8, vals=1:ncell(r))
r <- r * 1:8
time(r) <- as.Date(c('2020-01-01', '2020-01-02', '2020-02-01', '2020-02-02', '2021-01-01', '2021-01-02',  '2021-02-01', '2021-02-02'))

得到每年/每月的总和，然后求平均值，得到跨年的月平均值。

ym <- tapp(s, "yearmonths", sum)
m <- tapp(ym, "months", mean)

m
class       : SpatRaster 
dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source(s)   : memory
names       : m_1,  m_2 
min values  :   7,   11 
max values  : 700, 1100 
time (mnts) : Jan to Feb 

time(m, "months")
#[1] 1 2

Coy的解决方案运行良好（尽管他对问题的解释不同），但由于已经内置了日期操作，因此不必要地复杂。您可以使用以下公式获得他的结果：

y <- tapp(ym, "years", mean)
y
#class       : SpatRaster 
#dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory
#names       : y_2020, y_2021 
#min values  :      5,     13 
#max values  :    500,   1300 
#time (years): 2020 to 2021

如果你用rast(filename.nc)读取你的（netCDF）文件，你可能会得到自动设置的时间戳，这样你就不需要像我在示例数据中那样设置它们。

#nms <- names(x)
nms <- c("X1971.01.01",  "X1971.01.02")
as.Date(nms, "X%Y.%m.%d")
#[1] "1971-01-01" "1971-01-02"