Повысить производительность / скорость

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

Первый бит может быть просто:

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)

Вы также можете ускорить вторую часть, как показано в следующих примерах, в зависимости от того, сколько у вас ОЗУ.

Читать строку за строкой:

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 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 с против 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 перед выполнением цикла растров.Если вы соединяете 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 раз быстрее.stack поэтому я не могу сказать вам, как применить это к вашему коду.

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))
    }
  }
}