Augmenter la performance/vitesse

Question

J'ai besoin de prendre des données à partir de 1303 rasters (chaque raster a des données pendant 1 mois) et de faire une série chronologique pour chaque cellule de la grille dans les rasters. À la fin, je vais rejoindre toutes les séries chronologiques dans un fichier massif (zoo).

J'ai le code qui peut le faire (j'ai essayé sur une petite partie de l'ensemble de données et cela a fonctionné) mais il semble prendre pour toujours juste pour empiler le raster (plus de 2 heures maintenant et comptant toujours) ce n'est pas la partie la plus lente, qui fera la série chronologique. Donc voici mon code, si quelqu'un connaît un moyen plus rapide d'empiler des rasters et/ou de créer des séries temporelles (peut-être sans la double boucle?) Merci de nous aider ...

Je ne connais pas d'autre langage de programmation mais Serait-ce trop demander à 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

Je reposter mon commentaire ici et donner un meilleur exemple:

L'idée générale: allouer l'espace pour s avant la « raster' boucle est exécutée. Si vous concaténez s et r vers un nouvel objet à l'intérieur de la boucle, R doit allouer de la nouvelle mémoire pour chaque itération. C'est vraiment lent, surtout si s est grand.

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 

Comme vous pouvez le voir, la pré-allocation est environ 10 fois plus rapide.

Malheureusement, je ne connais pas raster et stack donc je ne peux pas vous dire comment l'appliquer à votre code.

Answer 2

Quelque chose comme cela devrait fonctionner (si vous avez assez de mémoire):

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

Il charge tous raster objets dans une grande liste, puis une fois les appels stack.

Voici un exemple des gains de vitesse: 1,25 sec contre 8 secondes pour 60 fichiers - mais l'ancien code est quadratique dans le temps pour que les gains sont beaucoup plus élevés pour plus de fichiers ...

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

Le premier bit pourrait être simplement:

s <- stack(files) 

la raison pour laquelle la création d'une pile est un peu lent est que chaque fichier doit être ouvert et vérifié pour voir si elle a la même nrow, etc. ncol que les autres fichiers. Si vous êtes absolument certain qui est le cas, vous pouvez utiliser un raccourci comme celui-ci (généralement pas recommandé)

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) 

Vous pouvez accélérer probablement la deuxième partie comme dans les exemples ci-dessous, selon la quantité de RAM vous avez.

ligne par ligne Lire:

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

plus extrême, lire toutes les valeurs en une seule fois:

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

j'ai essayé une autre façon de traiter avec de nombreux fichiers. J'ai d'abord combiné le raster de séries temporelles en un fichier au format NetCDF, en utilisant write.Raster (x, format = "CDF", .. puis juste lire un fichier pour chaque année, cette fois j'ai utilisé la brique (netcdffile, varname = ''), la lecture économise beaucoup. Cependant, je dois enregistrer la valeur de chaque cellule pour toutes les années selon un format prédéfini, dans lequel j'utilise write.fwf (x = v, ..., append = TRUE) mais cela prend beaucoup de temps pour près de 500 000 points. Quelqu'un a-t-il les mêmes expériences et la même aide pour accélérer ce processus? Voici mon code pour extraire toute la valeur pour chaque point:

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