0

I just picked up R this summer.

I created a nested for loop so that nls() will be performed in 24 raster layers and a new file will be created to save the data.

In this example, I have a raster stack of 48 layers worth 2 years of data. My for loop will perform the calculation in each year which consists of 24 layers, write a new ENVI file and proceed to the next year.

I managed to create a files for i = 1 and i = 2. However, R does not stop after creating new ENVI files. It is loading for forever and the temporary files in my C drive is updating from time to time (C:\Users\CandiceYingWoei\AppData\Local\Temp\Rtmp is empty, but the date modified keeps updating).

I tried to use break to stop the for loop. I think it is not necessary because i will never go over nyears since I already set i in (1:nyears). if (i > nyears) { break }

How do I stop R after creating the files?

library(sp)
library(raster)
library(ncdf4)
library(rgdal)

setwd("C:/")
new <- stack("82_83_test.envi")
nlayers <- nlayers(new)
nyears <- nlayers(new)/24

for (i in 1:nyears) {
  for (j in (1+(24*(i-1))):(24*i)) {
     fun<- function(x) { 
        # a function that perform nls through each pixel in all raster 
        # layers in each year and return several parameters (SOS, EOS, etc).
        return (c(SOS,EOS,LOS,SPUDOY,P_TAmplitude))
        }
     equation<-calc(new[[(1+(24*(i-1))):(24*i)]],fun,forceapply=TRUE)
     date<- 1981+i
     name<-paste("Output_",date,".envi",sep="")
     writeRaster(equation,filename=name,format="ENVI",overwrite=T)  
     # this is where the break statement is located
     if (i > nyears) { break }
  }
}
0

If you have to break a loop under this kind of circumstance, it should be a signal that you need to try something else. In this case, using lists and lapply would be a better idea. Try this approach instead, which will break your raster stack into a list of sub-stacks, one for each year, which you can then process in turn:

# just playing with some random files on this machine
setwd('C:/DATA')
x <- list.files(getwd(), pattern = '^MV_NSG.*\\.tif$', full.names = TRUE)

mvs <- stack(x)
mvlist <- unstack(mvs) # now a list of rasters (I have 6)
grpsize <- 3 # desired layers per stack - 12 for you

# use names to group - here, names(mvlist) becomes c('1','1','1','2','2','2')
# length.out will help if you have, say an extra 3 months of data in the original stack
names(mvlist) <- rep(1:(ceiling(length(mvlist) / grpsize )), 
                     each = grpsize, length.out = length(mvlist))

# make list of rasters into a list of stacks
# this gives me a list of 2 stacks with three rasters each 
stacks <- lapply(unique(names(mvlist)), function(y) {
  b <- mvlist[names(mvlist) == y]
  stack(b)
})

# now process the list - w is each stack in turn
stacksout <- lapply(stacks, function(w) {

  g <- calc(w, fun, forceapply = TRUE, etc)

})

consider using clusterR() to speed up your calc function too.

  • Thank you very much for your help @obrl_soil! My loop actually ended but it took a few hours for R to stop loading, after the raster files are written. I will try your method to see what works best. I do consider lapply as an option, however, it uses a list, so I need to convert the raster layers into a list, as what you do here, before performing the calculation, which I think could be less efficient. – Candice Jun 14 '17 at 20:57
0

You are start by calling year 1, then the j loop runs 24 times. Then repeats for year 2.

you could put the required layer numbers in a list, then loop through list elements:

nlayers <- nlayers(st)

l <- list(1:nlayers(new)/2,(nlayers(new)+1):14)

for (i in 1:length(l)){

    fun<- function(x) { 
        # a function that perform nls through each pixel in all raster 
        # layers in each year and return several parameters (SOS, EOS, etc).
        return (c(SOS,EOS,LOS,SPUDOY,P_TAmplitude))
        }
     equation<-calc(new[[l[[i]]]],fun,forceapply=TRUE)
     date<- 1981+i
     name<-paste("Output_",date,".envi",sep="")
     writeRaster(equation,filename=name,format="ENVI",overwrite=T)  

}
  • Thanks for the suggestion @Sam. I don't really understand the second line of your command. So you assign a list to variable l which in your case is (1:24, 49:14), assuming nlayers(new) = 48. But what is the purpose of this line? Also, in the calc() function you are calculating new[[l[[i]]]], which is just one raster layer. However, I would like to calculate it for the whole stack 24 layers for a year and do the same thing for the next year and so forth. – Candice Jun 14 '17 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.