使用R计算NetCDF数据的月平均气候统计量

Question

我正在处理为期31年的NOAA AVHRR每日海表面温度（SST）数据。该数据以NetCDF格式呈现，其维数为28（经度）x 40（纬度）x 11686（天）。我需要计算月度气候平均值（例如31年所有1月份的平均值，以此类推）。使用ncdf4和chron库，我能够将其转化成数组形式。

ncin <- nc_open('sstfile.nc')
sst_array <- ncvar_get(ncin, 'sst')

由于时间变量与SST数据是分开的，所以我必须使用它循环遍历数组。

is.leapyear <- function(year){
return(((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0))
}

dateseq <- seq(as.Date("1987-01-01"), as.Date("2018-12-31"), by=1)

使用光栅库将其转换为栅格图像，然后进行计算。

for ( i in seq(11686)) {
dtft <- strsplit(as.character(as.Date(dateseq[i])), split = '-')
y <-  as.integer(dtft[[1]][1])
m <-  as.integer(dtft[[1]][2])
d <-  as.integer(dtft[[1]][3])
while (m == 1){
assign(paste0('r',y,'.',d), raster(matrix(sst_array[1:27, 1:38, i], 
nrow = 27, ncol = 38)))
m = m + 1
}
if (is.leapyear(y) == TRUE) (i = i + 366)
else (i = i + 365)
}

问题在于它创建了太多栅格数据，首先计算月平均值，然后再计算年度平均值。

r87jan <- stack(mget(paste0('r1987.',1:31)))
r87janmean <- calc(r87jan, mean)

有没有可以在不生成过多栅格的情况下计算这段时间的函数/方法，使得计算结果可以保留为数组或矩阵？或者上述代码能否改进，一次性计算所有年份的月平均值？

- Satish_P

Answer 1

如果您已经安装了cdo（气候数据操作器），并且想要使用非R语言来回答此问题，您可以在Linux命令行上执行以下操作：

cdo ymonmean sstfile.nc sst_climate.nc 

文件sst_climate.nc将包含12个时间步长，分别为所有1月份、2月份等的平均值...

您可以在Ubuntu/Mint等操作系统中轻松安装cdo：

sudo apt-get install cdo 

现在你可以在Windows 10上轻松安装Ubuntu，以便轻松访问这些有用的工具。文档在这里https://code.mpimet.mpg.de/projects/cdo/

- Adrian Tompkins

Answer 2

您没有提供数据，但我认为您可以像这样做：

library(raster)
nc <- brick('sstfile.nc')

dates <- getZ(nc)
months <-  as.integer(format(dates, "%m"))

s <- stackApply(nc, months, fun=mean)

- Robert Hijmans