2

1303 個のラスター (各ラスターには 1 か月分のデータがあります) からデータを取得し、ラスター内の各グリッド セルの時系列を作成する必要があります。最後に、すべての時系列を 1 つの大規模な (動物園) ファイルに結合します。

私はそれを実行できるコードを持っています（データセットのごく一部を試してみましたが、うまくいきました）が、ラスターをスタックするだけで永遠にかかるようです（現在2時間以上、まだカウント中です）。より遅い部分で、時系列を実行します。だからここに私のコードがあります.誰かがラスターを積み重ねたり、時系列を作成したりするためのより速い方法を知っているなら(おそらく二重ループなしで?)助けてください...

私は他のプログラミング言語を知りませんが、これはRに尋ねるには多すぎますか? 

files <- list.files(pattern=".asc") 
pat <- "^.*pet_([0-9]{1,})_([0-9]{1,}).asc$"
ord_files <- as.Date(gsub(pat, sprintf("%s-%s-01", "\\1", "\\2"), files))
files<-files[order(ord_files)]


#using "raster" package to import data 
s<- raster(files[1])
pet<-vector()
for (i in 2:length(files))
{
r<- raster(files[i])
s <- stack(s, r)
}

#creating a data vector
beginning = as.Date("1901-01-01")
full <- seq(beginning, by='1 month', length=length(files))
dat<-as.yearmon(full)

#building the time series
for (lat in 1:360)
for (long in 1:720)
{
pet<-as.vector(s[lat,long])
x <- xts(pet, dat)
write.zoo(x,file=paste("P:/WRSRL/Users1/ncgk/IBERIA/cru_pet/zoo/","lat",lat,"long",long,".csv", sep="")  , sep=",")
}
4

4 に答える 4

2

最初のビットは単純に次のようになります。

s <- stack(files)

スタックの作成がやや遅い理由は、各ファイルを開いて、他のファイルと同じ nrow、ncol などがあるかどうかを確認する必要があるためです。その場合は、このようなショートカットを使用できます (一般的にはお勧めしません)。

quickStack <- function(f) {
r <- raster(f[1])
ln <- extension(basename(f), '')
s <- stack(r)
s@layers <- sapply(1:length(f), function(x){ r@file@name = f[x]; r@layernames=ln[x]; r@data@haveminmax=FALSE ; r })
s@layernames <- ln
s
}

quickStack(files)

RAM の容量によっては、次の例のように 2 番目の部分を高速化することもできます。

行ごとに読む:

for (lat in 1:360) {
pet <- getValues(s, lat, 1)
for (long in 1:720) {
    x <- xts(pet[long,], dat)
    write.zoo(x,file=paste("P:/WRSRL/Users1/ncgk/IBERIA/cru_pet/zoo/","lat",lat,"long",long,".csv", sep="")  , sep=",")
}
}

より極端な場合は、すべての値を一度に読み取ります。

 pet <- getValues(s)
 for (lat in 1:360) {
for (long in 1:720) {
    cell <- (lat-1) * 720 + long
    x <- xts(pet[cell,], dat)
    write.zoo(x,file=paste("P:/WRSRL/Users1/ncgk/IBERIA/cru_pet/zoo/","lat",lat,"long",long,".csv", sep="")  , sep=",")
}
}
于 2012-03-30T04:37:26.813 に答える
1

このようなものが動作するはずです (十分なメモリがある場合):

#using "raster" package to import data 
rlist <- lapply(files, raster)
s <- do.call(stack, rlist)
rlist <- NULL # to allow freeing of memory

すべてのオブジェクトを大きなリストにロードしてから、1 回raster呼び出しますstack。

速度向上の例を次に示します: 60 ファイルで 1.25 秒対 8 秒 - しかし、古いコードは時間的に 2 次であるため、より多くのファイルで速度向上がはるかに高くなります...

library(raster)
f <- system.file("external/test.grd", package="raster")
files <- rep(f, 60)

system.time({
 rlist <- lapply(files, raster)
 s <- do.call(stack, rlist)
 rlist <- NULL # to allow freeing of memory
}) # 1.25 secs

system.time({
 s<- raster(files[1])
 for (i in 2:length(files)) {
  r<- raster(files[i])
  s <- stack(s, r)
 }
}) # 8 secs
于 2012-03-29T18:32:14.033 に答える
0

多数のファイルを処理する別の方法を試しました。最初に、write.Raster(x,format="CDF",..) を使用して、時系列ラスターを NetCDF 形式の 1 つのファイルに結合し、次に、毎年 1 つのファイルを読み取るだけで、今回は brick(netcdffile,varname を使用) ='') 読み値が大幅に節約されます。ただし、 write.fwf(x=v,...,append=TRUE) を使用する事前定義された形式に従って、すべての年の各セルの値を保存する必要がありますが、500,000 ポイント近くには長い時間がかかります。このプロセスをスピードアップする方法について、同じ経験と助けを持っている人はいますか? 各ポイントのすべての値を抽出するためのコードは次のとおりです。

weather4Point <- function(startyear,endyear)  
{

  for (year in startyear:endyear)
  {
    #get the combined netCDF file

    tminfile <- paste("tmin","_",year,".nc",sep='')

    b_tmin <- brick(tminfile,varname='tmin')

    pptfile <- paste("ppt","_",year,".nc",sep='')

    b_ppt <- brick(pptfile,varname='ppt')

    tmaxfile <- paste("tmax","_",year,".nc",sep='')

    b_tmax <- brick(tmaxfile,varname='tmax')

    #Get the first year here!!!

    print(paste("processing year :",year,sep=''))

    for(l in 1:length(pl))
    {
      v <- NULL

      #generate file with the name convention with t_n(latitude)w(longitude).txt, 5 digits after point should be work

      filename <- paste("c:/PRISM/MD/N",round(coordinates(pl[l,])[2],5),"W",abs(round(coordinates(pl[l,])[1],5)),".wth",sep='')  

      print(paste("processing file :",filename,sep=''))            

      tmin <- as.numeric(round(extract(b_tmin,coordinates(pl[l,])),digits=1))

      tmax <- as.numeric(round(extract(b_tmax,coordinates(pl[l,])),digits=1))

      ppt <- as.numeric(round(extract(b_ppt,coordinates(pl[l,])),digits=2))

      v <- cbind(tmax,tmin,ppt)

      tablename <- c("tmin","tmax","ppt")

      v <- data.frame(v)   

      colnames(v) <- tablename

      v["default"] <- 0

      v["year"] <- year

      date <- seq(as.Date(paste(year,"/1/1",sep='')),as.Date(paste(year,"/12/31",sep='')),"days")

      month <- as.numeric(substr(date,6,7))

      day   <- as.numeric(substr(date,9,10))

      v["month"] <- month 

      v["day"]  <-  day

      v <- v[c("year","month","day","default","tmin","tmax","ppt")]

      #write into a file with format
      write.fwf(x=v,filename,append=TRUE,na="NA",rownames=FALSE,colnames=FALSE,width=c(6,3,3,5,5,5,6))
    }
  }
}
于 2013-10-31T15:17:50.857 に答える