2017-08-02 117 views
-1

我正在嘗試創建120個光柵文件的光柵堆棧。我在循環內生成這些文件,並將生成的柵格添加到rasterstack。代碼如下:堆棧函數僅堆棧最後一個光柵文件

library(raster) 
stack_P_95 <- stack() 
for (i in startyear:endyear) 
{ 
file <- paste(indir,"\\prec_",i,".nc",sep="") 
command <- paste("cdo timmin ",file," ",workdir,"min.nc",sep="") 
system(command) 
command <- paste("cdo timmax ",file," ",workdir,"max.nc",sep="") 
system(command) 
command <- paste("cdo timpctl,95 ",file," ",workdir,"min.nc 
",workdir,"max.nc ",workdir,"P95_",i,".nc",sep="") 
system(command) 
grid <- raster(paste(workdir,"P95_",i,".nc",sep="")) 
stack_P_95 <- stack(stack_P_95,grid) 
} 
crs(stack_P_95) <- "+proj=utm +zone=45 +ellps=WGS84 +datum=WGS84 +units=m 
+no_defs" 
writeRaster(stack_P_95,paste(outdir,model,"_P95.nc",sep=""),format="CDF",overwrite=T) 

但是,我創建只是有最後柵格重複n的stack_P_95(endyear -startyear)次。

什麼可能是這個問題的可能原因? 在此過程中沒有生成錯誤或警告。

我已經添加了鏈接與我輸入文件和最終產出: Input files

+0

你可以添加一個可重複的例子。例如'base'中沒有叫'addLayer'的函數。 – drmariod

+1

適用於我。你確定你的數據文件不完全相同嗎? – Spacedman

+0

下面的代碼用於創建當前目錄中包含年份數字網格的光柵文件測試集:'r = raster(); for(i in 2001:2004){r [] = i; writeRaster(r,paste0 「temp _」,i,「。nc」))}'#WARNING可以覆蓋你的文件#這樣做,然後運行你的代碼,我會得到一個4層堆棧,其中包含四層數值。即它適用於我。 – Spacedman

回答

0

沒有什麼錯與您的代碼。

創建一些示例柵格:

> library(raster) 
Loading required package: sp 
> for(i in 2001:2005){ r = raster(matrix(i,10,10)); writeRaster(r,paste0("temp_",i,".nc"))} 
Loading required namespace: ncdf4 

這給了我:

$ ls 
temp_2001.nc temp_2002.nc temp_2003.nc temp_2004.nc temp_2005.nc 

所有這些都是不同的:

$ gdalinfo temp_2004.nc | grep max 
    layer#max=2004 
    max=2004 
$ gdalinfo temp_2001.nc | grep max 
    layer#max=2001 
    max=2001 

現在在清潔R對話我跑你的代碼:

library(raster) 
workdir="./" 
stack_1 <- stack() 
startyear = 2001 
endyear = 2004 
for (i in startyear:endyear) 
{ 
tempnc <- raster(paste(workdir,"temp_",i,".nc",sep="")) 
stack_1 <- stack(stack_1,tempnc) 
} 

我也得到了四種不同層的堆疊:

> range(stack_1[[1]][]) 
[1] 2001 2001 
> range(stack_1[[2]][]) 
[1] 2002 2002 
> range(stack_1[[3]][]) 
[1] 2003 2003 

不像你的代碼,我的代碼是完全可重複的。如果你在一個乾淨的R環境中運行我的代碼,並且所有圖層都是「2004」,那麼你的系統出現了一些可怕的錯誤。如果圖層與我的圖層相同,那麼您的文件或代碼不會顯示給我們。

+0

我已將我的整個代碼以及我的輸入文件和輸出文件添加到問題中。我正在使用一個操作的輸出,並將它讀回堆棧,這在我的情況下不起作用。 – Anonymousaurav

0

試試這個代碼

r <- list() # create list objects 

for (i in startyear:endyear) 
{ 
r[i] <- raster(paste(workdir,"temp_",i,".nc",sep="")) 
} 

rs <- stack(r) # create raster stack from list object 
+0

如果我第一次生成輸出,然後按照您的建議列出它們,堆疊工作。但我不明白爲什麼堆棧在循環中不起作用。 – Anonymousaurav