基于堆叠的NC日降水数据，如何计算各经纬度栅格月均值？

我明白你的问题了——你需要的是每个经纬度单元格在对应年份1月的逐日降水均值，而不是整个栅格层的单一全局均值。CellStats()函数的作用是计算整个栅格的统计量（比如所有单元格的平均值），所以它返回单一数值是正常的。你需要的是对每个单元格跨时间维度（该年份1月的所有日期）计算均值，下面分两种常用的R包场景给出具体实现：

场景1：使用raster包（适配你参考的合并教程） #

假设你已经将每个年份1月的逐日数据堆叠成了单独的RasterStack（比如jan_1981对应1981年1月的31个逐日栅格），或者所有年份1月的逐日数据合并成了一个大的RasterStack：

情况A：每个年份1月是单独的栅格栈 #

对单个年份的1月栈，使用calc()函数对每个单元格应用均值计算：

library(raster)

# 以1981年1月为例
jan_1981_stack <- stack("1981_jan_daily.nc") # 或你已读入的栈对象
jan_1981_mean <- calc(jan_1981_stack, fun = mean, na.rm = TRUE)

# 保存结果（可选）
writeRaster(jan_1981_mean, "1981_jan_mean.nc", format = "CDF")

重复这个步骤处理所有年份的1月栈即可。

情况B：所有年份1月的逐日数据在同一个大栈里 #

如果你的大栈图层名称包含年份信息（比如"1981-01-01"、"1981-01-02"...），可以用stackApply()按年份分组批量计算：

all_jan_daily <- stack("all_jan_daily_merged.nc")

# 从图层名称提取年份
year_indices <- substr(names(all_jan_daily), 1, 4)

# 按年份分组，对每组（每个年份1月的所有日期）计算单元格均值
all_jan_monthly <- stackApply(all_jan_daily, indices = year_indices, fun = mean, na.rm = TRUE)

# 保存所有年份的1月均值栈
writeRaster(all_jan_monthly, "all_jan_monthly_means.nc", format = "CDF")

场景2：使用terra包（推荐，raster包已停止维护） #

terra是raster的替代包，效率更高，语法更直观：

情况A：单个年份1月的逐日数据 #

library(terra)

# 读取单个年份1月的逐日数据
jan_1981 <- rast("1981_jan_daily.nc")

# 对每个单元格计算1月均值
jan_1981_mean <- app(jan_1981, mean, na.rm = TRUE)

# 保存结果
writeCDF(jan_1981_mean, "1981_jan_mean.nc", overwrite = TRUE)

情况B：所有年份1月的逐日数据在一个对象里 #

如果你的栅格对象带有时间属性，可以直接用tapp()按年份分组：

all_jan_daily <- rast("all_jan_daily_merged.nc")

# 从时间属性提取年份
year_indices <- format(time(all_jan_daily), "%Y")

# 按年份分组计算月均值
all_jan_monthly <- tapp(all_jan_daily, index = year_indices, fun = mean, na.rm = TRUE)

# 保存结果
writeCDF(all_jan_monthly, "all_jan_monthly_means.nc", overwrite = TRUE)

关键区别再强调一下：CellStats()是全局统计（整个栅格的单一值），而calc()/app()/stackApply()/tapp()是逐单元格统计（每个位置保留一个均值，生成新的栅格）。

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