Leistung/Geschwindigkeit erhöhen

Question

Ich muss Daten aus 1303 Rastern nehmen (jedes Raster hat Daten für 1 Monat) und eine Zeitreihe für jede Rasterzelle in den Rastern erstellen. Am Ende werde ich alle Zeitreihen zu einer riesigen (Zoo-) Datei zusammenfügen.

Ich habe den Code, der es tun kann (Ich versuchte einen kleinen Teil des Datasets und es funktionierte), aber es scheint für immer nur das Raster zu stapeln (mehr als 2 Stunden jetzt und immer noch gezählt) und das ist nicht der langsamere Teil, der die Zeitreihe machen wird. Also hier ist mein Code, wenn jemand einen schnelleren Weg zum Stapeln von Rastern kennt und/oder um die Zeitreihe zu erstellen (vielleicht ohne Doppelschleife?) Bitte helfen ...

Ich kenne keine andere Programmiersprache aber Wäre das zu viel von R zu fragen?

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

Ich werde meinen Kommentar hier umbuchen und ein besseres Beispiel geben:

Die allgemeine Idee: zuteilen den Raum für s, bevor die ‚raster'-Schleife ausgeführt wird. Wenn Sie s und r mit einem neuen Objekt s innerhalb der Schleife verketten, muss R für jede Iteration neuen Speicher für s zuweisen. Das ist sehr langsam, besonders wenn s groß ist.

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 

Wie Sie sehen können, ist die Vorbelegung etwa 10 mal schneller.

Leider bin ich nicht vertraut mit raster und stack so kann ich Ihnen nicht sagen, wie Sie dies auf Ihren Code anwenden.

Answer 2

So etwas sollte funktionieren (wenn Sie über genügend Speicher):

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

Es lädt alle raster Objekte in eine große Liste und ruft dann stack einmal.

Hier ist ein Beispiel für die Geschwindigkeit gewinnt: 1,25 sec vs 8 Sekunden für 60 Dateien - aber Ihr alter Code ist in der Zeit quadratisch, so dass die Gewinne für mehr Dateien viel höher sind ...

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

Die erste Bit könnte einfach sein:

s <- stack(files) 

Der Grund, warum die einen Stapel zu schaffen etwas langsam ist, dass jede Datei geöffnet werden muss und geprüft, um zu sehen, ob es die gleiche nrow hat, ncol usw. wie die anderen Dateien. Wenn Sie absolut sicher das der Fall ist, können Sie eine Verknüpfung wie folgt verwenden (im Allgemeinen nicht empfohlen)

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) 

Sie können sich wahrscheinlich auch den zweiten Teil, wie in den folgenden Beispielen beschleunigen, je nachdem, wie viel RAM hast du.

liest Zeile für Zeile:

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

extreme, lesen Sie alle Werte in einem Rutsch:

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

habe ich versucht, eine andere Art und Weise mit einer Vielzahl von Dateien zu tun hat. Zuerst habe ich das Zeitreihen-Raster in eine Datei im NetCDF-Format, , kombiniert. Mit write.Raster (x, format = "CDF", ..) und dann nur eine Datei für jedes Jahr zu lesen, dieses Mal habe ich Ziegel (netcdffile, varname = '') verwendet es spart viel. Allerdings muss ich den Wert jeder Zelle für alle Jahre nach einem vordefinierten Format speichern, in dem ich write.fwf (x = v, ..., append = TRUE) , aber es dauert eine lange Zeit für fast 500.000 Punkte. Hat jemand die gleichen Erfahrungen und hilft, diesen Prozess zu beschleunigen? Hier ist mein Code alle den Wert für jeden Punkt zum Extrahieren:

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