提高性能/速度

Question

我需要从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=",") 
} 

Answer 1

我将在这里重新发布我的意见，并给一个更好的例子：

总体思路：在执行“raster'环之前分配对于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，所以我不能告诉你如何将它应用到你的代码中。

Answer 2

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

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

它加载所有raster物体进入大名单，然后调用stack一次。

这里的速度上涨的一个例子：1.25秒VS 8秒60页的文件 - 但你的旧代码是二次在时间上的收益更多的文件要高得多......

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 

Answer 3

第一位可能仅仅是：

s <- stack(files) 

之所以创建堆栈是有点慢是每个文件需要被打开和检查，看它是否具有相同的nrow，NcoI位等的其他文件。如果你是绝对肯定是这样的话，你可以使用这样的快捷方式（一般不推荐）

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=",") 
} 
} 

Answer 4

我尝试另一种方式来处理大量文件。 首先，我将时间序列栅格合并到NetCDF格式的一个文件中， 使用write.Raster（x，format =“CDF”，..） 然后每年只读一个文件，这次我使用了砖块（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)) 
    } 
    } 
}