2

したがって、4D 配列 (x、y、z、t) にある .nc ファイルの変数がいくつかあります。問題は、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)  
        }  
    }  
}

これをより速く行う方法はありますか？または、これを簡単にするためにデータを再グリッド化することを検討するよう誰かに言われましたが、彼が何を意味しているのかよくわかりません。誰かが私を助けることができますか？ありがとう。

4

2 に答える 2

3

あなたの 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は、値のベクトルも取ることができます。ここではapprox、次元 1、2、および 4 に関数を適用します (変更するのは 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

時間の増加は比較的小さいですが、おそらくそれだけの価値があります。

于 2014-10-20T08:06:00.953 に答える