使用 R 估计 NetCDF 数据的每月气候
我正在研究 31 年的 NOAA AVHRR 每日海面温度 (SST) 数据.数据采用 NetCDF 格式,尺寸为 28(lon)x 40(lat)x 11686(天).我应该计算每月的气候平均值(例如,所有 31 年一月的平均值等等).使用 ncdf4 和 chron 库,我能够以数组形式获取它.
I am working on NOAA AVHRR daily Sea Surface Temperature (SST) data of 31 years. The data is in NetCDF format with dimensions as 28 (lon) x 40 (lat) x 11686(days). I am supposed to compute monthly climatological mean (e.g. mean of all Januaries of 31 years and so on). Using ncdf4 and chron libraries I was able to get it in array form.
ncin <- nc_open('sstfile.nc')
sst_array <- ncvar_get(ncin, 'sst')
由于时间变量与 SST 数据是分开的,我不得不使用它在数组上循环.
Since time variable is separate from SST data, I had to use it loop over the array.
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)
使用栅格库,我将转换为栅格,然后进行计算.
Using raster library I am converting to rasters and then doing computation.
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)
}
问题是它创建了太多的栅格,首先计算月平均值,然后每年计算.
Problem is it's creating far too many rasters and first computing monthly mean and then yearly.
r87jan <- stack(mget(paste0('r1987.',1:31)))
r87janmean <- calc(r87jan, mean)
是否有任何函数/方法可以在这段时间内进行计算而无需制作这么多栅格并且计算可以保留为数组或矩阵?或者可以改进上面的代码以一次计算所有年份的月平均值?
Is there any function/method which can compute over this time duration without making so many rasters and computation can remain as array or a matrix? Or can the above code can be improved to compute the monthly mean for all years at once?
你没有提供你的数据,但我认为你可以这样做:
You do not provide your data, but I think you can do something like this:
library(raster)
nc <- brick('sstfile.nc')
dates <- getZ(nc)
months <- as.integer(format(dates, "%m"))
s <- stackApply(nc, months, fun=mean)