R中的regrid netcdf数据用于插值答案

【问题标题】：regrid netcdf data in R for interpolationR中的regrid netcdf数据用于插值
【发布时间】：2014-10-20 15:50:58
【问题描述】：

所以，我有一些来自 .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)  
        }  
    }  
}

有没有办法更快地做到这一点？或者，有人告诉我研究重新划分数据以使这更容易，但我不太确定他的意思。有人能帮我吗？谢谢。

【问题讨论】：

您需要处理所有关卡吗？因为在我看来，你想要达到的目标没有意义。
是的，我需要所有关卡。简而言之，我正在做的是跟踪云和云顶，并保存在 z 中跨越 4 公里的数据矩阵。由于 z 坐标在地面附近更紧密地排列，这将导致矩阵大小不同（对于靠近地面的云，z 可能需要 60 个点，但对于 z 较高的云，只需要 40 个点）。

标签： r interpolation netcdf r-grid cdo-climate

【解决方案1】：

由于我没有你的 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函数（因为它是我们正在修改的第 3 个维度）：

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

时间增益相对较小，但可能值得。

【讨论】：

【解决方案2】：

这不是R中的答案，只是说这个任务可以用CDO从命令行快速完成

 cdo intlevel,`seq -s "," 50 50 22000` in.nc out.nc

seq 命令以 50m 的间隔生成一个从 50 到 22000 的逗号分隔列表。

【讨论】：