R语言计算逐日栅格周均值时的闰年天数适配问题求解

项目背景 #

我正在开展时序温度数据追踪相关项目，已生成覆盖40年时间范围的逐日栅格地图，后续计划计算逐周平均结果并导出为.png格式文件，最终合成GIF实现各年度数据的动态可视化。当前编写的代码在平年场景下可正常运行，但处理闰年数据时会出现运行故障。

原有实现代码 #

dseq <- seq(from = as.Date("01-01-1993", format = "%d-%m-%Y"), to = as.Date("31-12-1993", 
    format = "%d-%m-%Y"), by = 1)
hdu.fname <- paste("chdu", format(dseq, format = "%Y%m%d"), ".tif", sep = "")
img.fname <- paste("chdu_week_", 1:52, ".png", sep = "")
tif.fname <- paste("chdu_week_", 1:52, ".tif", sep = "")

wid <- c(rep(1:52, each = 7), 52)
ofile <- data.frame(wid, hdu = hdu.fname)

for(i in 1:52){
  id <- ofile$wid == i
  tofile <- ofile[id,]
  
  tStack <- stack()
  for(j in 1:length(tofile$hdu)){
     traster <- raster(as.character(tofile$hdu[j]))
     tStack <- stack(tStack, traster)
  }

  tchdu.r <- calc(tStack, mean)

  writeRaster(x = tchdu.r, filename = as.character(tif.fname[i]), overwrite = TRUE , format = "GTiff")
  
  breaks <- seq(from = 0, to = 800, by = 100)
  cols <- terrain.colors(n = length(breaks) - 1, alpha = 1)
  sdate <- as.character(tofile[1,2])
  
  sdate <- substr(sdate, start = 5, stop = 12)
  sdate <- as.Date(sdate, format = "%Y%m%d")
  sdate <- format(sdate, format = "%d-%b-%Y")

  windows(); plot(x = aupoly.ext[1:2], y = aupoly.ext[3:4], type = "n", xlab = "Longitude", ylab = "Latitude")
  image(tchdu.r, col = cols, breaks = breaks, add = TRUE, zlim = c(0,1000))
  plot(auadm0ll.sf, add = TRUE, colour="transparent", border="#696969")
  contour(tchdu.r, levels = 130, lty = 1, add = TRUE, lwd=1.5, col="purple", drawlabels=TRUE)
  text(x = aupoly.ext[1] + 4, y = aupoly.ext[3] + 2, labels = sdate)
  metre(xl = aupoly.ext[1], yb = aupoly.ext[3] + 24, xr = aupoly.ext[1] + 1, yt = aupoly.ext[3] + 33, lab = breaks, cols = cols, shift = 0, cex = 0.80)
  savePlot(filename = img.fname[i], type = c("png"), device = dev.cur())
  dev.off()
  
  cat(i, "\n"); flush.console()
}

故障表现 #

当将处理年份修改为闰年（例如1992年）时，dseq向量的长度为366，但wid向量的长度仅为365，二者行数不匹配导致无法创建ofile数据框，后续循环代码无法正常执行，触发报错如下：

Error in data.frame(wid, hdu = hdu.fname) :
arguments imply differing number of rows: 365, 366

故障原因 #

原有代码中wid的生成逻辑是硬编码实现：c(rep(1:52, each = 7), 52)，生成的序列固定长度为365，仅适配平年天数。遇到366天的闰年时，硬编码生成的周ID序列比实际日期序列少1个元素，直接导致数据框创建失败。

修复方法 #

• 移除硬编码的周ID生成逻辑，根据dseq的实际长度动态计算周编号，自动适配平年、闰年的天数差异
• 将原有wid <- c(rep(1:52, each = 7), 52)替换为以下代码：

# 按日期序列实际长度动态生成周ID
day_index <- seq_along(dseq)
wid <- ceiling(day_index / 7)
# 将年末不足7天的剩余日期统一归入第52周
wid[wid > 52] <- 52

替换后wid长度会和dseq、hdu.fname完全一致，不会再出现行数不匹配的报错，平年、闰年场景均可正常运行。

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