2

我需要从 1303 个栅格中获取数据（每个栅格都有 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。

逐行读取：

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

我将在这里重新发表我的评论并举一个更好的例子：

总体思路：在执行“光栅”循环之前为 s 分配空间。如果将 s 和 r 连接到循环内的新对象 s，则 R 必须为每次迭代为 s 分配新内存。这真的很慢，特别是如果 s 很大。

s <- c()
system.time(for(i in 1:1000){ s <- c(s, rnorm(100))})
# user  system elapsed 
# 0.584   0.244   0.885 

s <- rep(NA, 1000*100)
system.time(for(i in seq(1,1000*100,100)){ s[i:(i+99)] <- rnorm(100) })
# user  system elapsed 
# 0.052   0.000   0.050

如您所见，预分配快了大约 10 倍。

不幸的是我不熟悉，raster 所以stack我不能告诉你如何将它应用到你的代码中。

于 2012-03-29T17:54:55.217 回答
1

像这样的东西应该可以工作（如果你有足够的内存）：

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

它将所有对象加载raster到一个大列表中，然后调用stack一次。

这是速度增益的一个示例：60 个文件的速度为 1.25 秒与 8 秒 - 但是您的旧代码在时间上是二次的，因此对于更多文件来说，增益要高得多......

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 格式的文件中，然后每年只读取一个文件，这次我使用了 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 回答