我正在研究NOAA AVHRR每日海表温度(SST)数据31年。数据采用NetCDF格式,尺寸为28(长)x 40(纬度)x 11686(天)。我应该计算每月的气候平均值(例如31年所有Januaries的平均值等)。使用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)
是否有任何函数/方法可以在这段时间内计算而不需要制作如此多的栅格,并且计算可以保留为数组或矩阵?或者上面的代码可以更好地一次性计算所有年份?
您没有提供数据,但我认为您可以这样做:
library(raster)
nc <- brick('sstfile.nc')
dates <- getZ(nc)
months <- as.integer(format(dates, "%m"))
s <- stackApply(nc, months, fun=mean)
如果我能提供不在R中的答案,如果你有cdo(安装了气候数据操作员),你可以在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/