CODEWITHSUNDEEP
rspatialrasterr-raster
Candice
Candice

Reputation: 23

R: Read raster layers as one object to apply a function in each ENVI file

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!

Upvotes: 0

Views: 492

Answers (2)

Candice
Candice

Reputation: 23

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

Upvotes: 0

maRtin
maRtin

Reputation: 6516

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]]

Upvotes: 0

Related Questions