0

I have an ENVI file (82_83_test.envi) that contains biweekly raster layers from 1982 to 1983. That is 24 layers per year, 48 layers in total. I would like to create a for loop to apply a function to perform time-series analysis per year, i.e., R will run through 24 layers in a pixel and calculate 5 parameters with the function "fun" for all pixels for that year. Ultimately, I would like to have 5 plots (5 parameters) for each year, so a total of 10 plots for two years.

I tried working with 1 ENVI file with 2 years of data and 2 ENVI files with 1 year of data in each file. I used brickstack_to_raster_list() from the library spatial.tools to read the file, and I get 48 layers. However, I would like to get 2 chunks (1982 and 1983) which consist of 24 layers for each chunk, so that I can run the equation.

Maybe something like brickstack_to_raster_list() then merge the 1st layer to 24th into one, followed by the 25th to 48th layer into one?

new <- stack("82_83_test.envi")
new1<- brickstack_to_raster_list(new)

new1 returns 48 raster layers. For example,

new1
[[1]]
class       : RasterLayer 
band        : 1  (of  48  bands)
dimensions  : 151, 101, 15251  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -105.0833, -96.66667, 56.66667, 69.25  (xmin, xmax, ymin, 
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 
+towgs84=0,0,0 
data source : C:\*\82_83_test.envi 
names       : Band.1 
values      : -32768, 5038  (min, max)

The other approach is to concatenate multiple annual ENVI files into a list. 

new <- stack("1982_test.envi")
new1<- stack(new,new)
new2<- brickstack_to_raster_list(new1)

Both methods above yield the same result, although I am not certain about its efficiency. Because after getting this set up, I will be generating data from 1982 to 2015, and so efficiency matters a lot.

Below is the function that I would like to apply in the for loop.

# A is an unknown that will be the number of components in the list. 
for (i in length(A)) {
new1[new1<=-1000]<-0
Data_value<-new1/10000
# assign 0 to pixel value that is less than -1000 and divide by 10000 in order to use the equation
DOY<-(1:nlayers(new1)*15)
# so that the unit will be in days instead of the number of weeks.

fun<- function(x) { if (all(is.na(x[1]))) { return(rep(NA,5)) } else { 
fitForThisData  <-nls(x~ a+((b/(1+ exp(-c*(DOY-e))))- (g/(1+ exp(-d*(DOY-
f))))), alg="port",start=list(a=0.1,b=1,g=1,c=0.04,d=0.04,e=112,f=218),
lower=list(a=0,b=0.3,g=0.3,c=-1,d=-1,e=20,f=100),
upper=list(a=0.4,b=2,g=2,c=1,d=1,e=230,f=365),
control=nls.control(maxiter=2000, tol = 1e-15, minFactor = 1/1024, 
warnOnly=TRUE))
SOS<-(coef(fitForThisData)[6] -(4.562/(2*coef(fitForThisData)[4])))
EOS<-(coef(fitForThisData)[7] -(4.562/(2*coef(fitForThisData)[5])))
LOS<-(EOS-SOS)
SPUDOY<-(1.317*((-1/coef(fitForThisData)[4])+ coef(fitForThisData)[6]))
P_TAmplitude<-(SPUDOY-SOS)
return (c(SOS,EOS,LOS,SPUDOY,P_TAmplitude))
}
}
}

equation<-calc(Data_value,fun,forceapply=TRUE)
plot(equation)

I would truly appreciate your advice on how to do this. Thank you very much!

Candice
  • 31
  • 8

2 Answers2

0

After reading in your stack:

library(raster)    
new <- stack("82_83_test.envi")

simply split the stack into yearly substacks using basic indexing:

year1 <- new[[1:24]]
year2 <- new[[25:48]]
maRtin
  • 6,336
  • 11
  • 43
  • 66
  • Indexing sounds like a good way to do it. What about doing a nested for loop? Or is there another more efficient way of doing this? – Candice Jun 05 '17 at 20:43
  • Hi @maRtin, I tried indexing but I couldn't let writeRaster to write two separate files for me for different years. year1 <- new[[1:24]] year2 <- new[[25:48]] years<-c(year1,year2) for (i in 1:length(years)){ for (j in 1:nlayers(years[[i]])) { function (x) equation<-calc(Data_value,fun,forceapply=TRUE) date<- 1981+i writeRaster(equation,filename=paste("YW_Output_0606_",date,".envi",sep=""),format="ENVI",overwrite=T) } } The calc() was performed for two years, instead of each year. And the command writeRaster only writes one file, instead of two. Thoughts? Thank you. – Candice Jun 06 '17 at 19:49
  • @Candice Sorry, but it's really hard to replicate and understand what the problem is. A minimal reproducible example: https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example would be very helpful. Maybe you can eliminate code that works and construct a short example of the part thats not working and how the output should look like. I am sure that your question will attract more answers if you follow these guidelines. – maRtin Jun 07 '17 at 19:55
  • The example should include the follwoing: -) a minimal dataset, necessary to reproduce the error -) the minimal runnable code necessary to reproduce the error, which can be run on the given dataset. – maRtin Jun 07 '17 at 19:56
0

UPDATE: I was able to loop the function but my guess is that the calculation was done on all raster layers after comparing it with the truth value. However, two files of the same content with different file names are produced, because two files have the same summaries.

new <- stack("82_83_test.envi")
new[new<=-1000]<-0
Data_value<-new/10000
nlayers <- nlayers(new)
nyears <- nlayers(new)/24
DOY<-((1:nlayers(new))/nyears)*15
dummy<- FALSE

for (i in 1:nyears) {
  for (j in (1+24*(i-1)):(24*i)) { 
    fun<-function (x)
    equation<-calc(Data_value,fun,forceapply=TRUE)
    date<- 1981+i
    writeRaster(equation,filename=paste("Output",date,".envi",sep=""),
    format="ENVI",overwrite=T)
    if (j == nlayers){
    dummy<-TRUE
    break
    if (dummy) {break}      
}
}
Candice
  • 31
  • 8