如何对3D数组日数据计算月均值并存储为4D数组？报错求助

我有两年的每日地理空间数据，分别存储在两个NC文件中，每个文件包含365天的数据。需要先对数据执行自定义计算，再为每个网格单元计算每年12个月的月均值，但运行代码时出现下标越界错误，尝试rowMeans也未得到正确结果。

原始代码（出错版本） #

lons <- 584
lats <- 712
MyValues_ym <- array(0., c(lons, lats, 2, 12))
MyCalculation <- array(0., c(lons, lats, ts))  ## 初始化存储计算结果的数组
MyMonths = c('01','02','03','04','05','06','07','08','09','10','11','12')
ts <- dim(P)[[3]]  ## 匹配降水数据的时间维度（365天）
## 降水数据维度：[584L, 712L, 365L]

## 计算月均值
for (m in 1:12) {
  MyValues_ym[, , y - 1982, m] <- apply(MyCalculation, c(2), mean)
}

报错信息 #

Error in `[<-`(`*tmp*`, , , y - 1982, m, value = c(NA_real_, NA_real_,  : 
  subscript out of bounds

尝试用rowMeans的代码（未正确计算） #

# 循环计算月均值
for (m in 1:12) {     
  my_months_yearly <- format(daterange, '%m') == MyMonths[m]    
  MyValues_ym[,,y-1982,m] = rowMeans(Q0d[,,my_months_yearly], na.rm = TRUE, dims = 2) 
}

最小可复现示例（MRE） #

lons <- 10
lats <- 10

MyYearlyMonthlyData <- array(0., c(lons,lats,1,12)) ## 初始化月均值存储数组
MyCalculationYr1 <- array(runif(1:100), c(lons, lats, 365)) ## 模拟每日计算数据

## 设置日期范围和月份
start = '1985-01-01'
end = '1985-12-31'
daterange = seq(as.Date(start), as.Date(end), "days")
MyMonths = c('01','02','03','04','05','06','07','08','09','10','11','12')
y=1985

## 循环计算月均值
for (m in 1:12) { 
  blMM <- format(daterange,'%m') == MyMonths[m]
  MyYearlyMonthlyData[,,y-1984,m] = rowMeans(MyCalculationYr1[,,blMM],na.rm=TRUE,dims=2) 
}  

问题分析与修正方案 #

错误原因 #

1. 下标越界：MyValues_ym的第三个维度长度为2（对应两年数据），但y - 1982如果y取1985会得到3，超出了数组维度的有效索引（1-2），导致报错。
2. 均值计算维度错误：原始代码中apply(MyCalculation, c(2), mean)是对纬度维度取均值，而非针对每个网格单元（lon,lat）的时间维度取均值。
3. MRE数据生成错误：array(runif(1:100), c(lons, lats, 365))生成的随机数数量（100）远小于数组所需的36500个元素，会导致维度不匹配。

修正后的MRE代码 #

lons <- 10
lats <- 10

# 初始化月均值存储数组：[经度, 纬度, 年份, 月份]
MyYearlyMonthlyData <- array(0., c(lons, lats, 1, 12)) 
# 生成正确长度的模拟每日计算数据
MyCalculationYr1 <- array(runif(lons*lats*365), c(lons, lats, 365)) 

# 设置日期范围和月份
start <- '1985-01-01'
end <- '1985-12-31'
daterange <- seq(as.Date(start), as.Date(end), "days")
MyMonths <- c('01','02','03','04','05','06','07','08','09','10','11','12')
y <- 1985

# 循环计算每个月的均值
for (m in 1:12) { 
  # 筛选对应月份的日期索引
  blMM <- format(daterange,'%m') == MyMonths[m]
  # rowMeans的dims=2：将前2个维度（经度、纬度）视为单个网格单元，对时间维度取均值
  MyYearlyMonthlyData[,, y-1984, m] <- rowMeans(MyCalculationYr1[,, blMM], na.rm = TRUE, dims = 2) 
}  

针对两年数据的完整修正代码 #

假设已读取两年数据并完成自定义计算，得到MyCalculation1（1982年，维度[584,712,365]）和MyCalculation2（1983年，维度[584,712,365]）：

lons <- 584
lats <- 712
# 初始化存储数组：[经度, 纬度, 年份(2年), 月份(12)]
MyValues_ym <- array(0., c(lons, lats, 2, 12)) 
MyMonths <- c('01','02','03','04','05','06','07','08','09','10','11','12')

# 处理1982年数据（对应数组第3维度索引1）
daterange1 <- seq(as.Date('1982-01-01'), as.Date('1982-12-31'), "days")
for (m in 1:12) {
  blMM <- format(daterange1,'%m') == MyMonths[m]
  MyValues_ym[,, 1, m] <- rowMeans(MyCalculation1[,, blMM], na.rm = TRUE, dims = 2)
}

# 处理1983年数据（对应数组第3维度索引2）
daterange2 <- seq(as.Date('1983-01-01'), as.Date('1983-12-31'), "days")
for (m in 1:12) {
  blMM <- format(daterange2,'%m') == MyMonths[m]
  MyValues_ym[,, 2, m] <- rowMeans(MyCalculation2[,, blMM], na.rm = TRUE, dims = 2)
}

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