因此,我有一些来自 .nc 文件的变量,它们位于 4D 数组 (x,y,z,t) 中。问题是,z 坐标不像 x 和 y 坐标那样均匀分布,即 z 大约为 25 米、75m、125、175、...、500、600、700、...、20000、 21000, 22000。我试图对数据进行线性插值以获得整个 z 上均匀的 50m 间距。但是 R 中的 approx 函数运行得太慢(我认为数组太大):
library(ncdf)
x = get.var.ncdf(nc,'x'); y = get.var.ncdf(nc,'y'); z = get.var.ncdf(nc,'z')
t = get.var.ncdf(nc,'t') # time
qc1 = get.var.ncdf(nc,'qc',start=c(1,1,1,1),count=c(-1,-1,-1,-1))
zlin = seq(z[1],z[length(z)],50)
qc1_lin = array(0,c(length(x),length(y),length(zlin),length(t)))
for (i in 1:length(x)) {
for (j in 1:length(y)) {
for (k in 1:length(t)) {
qc1_lin[i,j,,k] = approx(z,qc1[i,j,,k],xout = zlin)
}
}
}
有没有办法可以更快地做到这一点?或者,有人告诉我研究重新网格化数据以使其更容易,但我不太确定他的意思。有人能帮我吗?谢谢。
最佳答案
由于我没有您的 ncdf 文件,因此我使用 NOAA 气温数据集作为示例:
library(ncdf)
url <- paste("ftp://ftp.cdc.noaa.gov/Datasets/ncep/air.",format(Sys.Date(),"%Y"),".nc",sep="")
download.file(url,destfile="air.nc")
nc <- open.ncdf("air.nc")
x <- get.var.ncdf(nc,'lon')
y <- get.var.ncdf(nc,'lat')
z <- get.var.ncdf(nc,'level')
t <- get.var.ncdf(nc,'time')
qc1 <- get.var.ncdf(nc,'air')
这里的z
值范围从1000到50,举一个简短的例子,让我们采用每100级间隔的规则网格(我还将限制数据集前20天的操作)以保持示例相对较小):
zlin <- seq(z[1],z[length(z)],-100)
使用你的方法:
qc1_lin <- array(0,dim=c(144,73,10,20))
system.time({
for (i in 1:length(x)) {
for (j in 1:length(y)) {
for (k in 1:20) {
# Don't forget that approx outputs a list
qc1_lin[i,j,,k] = approx(z,qc1[i,j,,k],xout = zlin)$y
}
}
}
})
user system elapsed
26.793 1.196 27.886
但是您可以使用 apply
执行相同的操作:参数 MARGIN
也可以采用值向量。这里我们要在维度 1、2 和 4 上应用approx
函数(因为我们正在修改的是第三个维度):
system.time({
qc1_lin2 <- apply(qc1[,,,1:20],c(1,2,4),function(X)approx(z,X,xout=zlin)$y)
})
user system elapsed
24.413 0.144 24.408
apply
不幸的是,输出新维度作为第一个维度,因此我们需要排列结果:
qc1_lin3 <- aperm(qc1_lin2, perm=c(2,3,1,4))
让我们检查一下结果是否相同:
all(qc1_lin3==qc1_lin)
[1] TRUE
时间增益相对较小,但可能是值得的。
关于R中的regrid netcdf数据用于插值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/25520755/