0

I am working with a time-series raster brick. The brick has 365 layers representing a value for each day of the year.

I want to create a new layer in which each cell holds the number of day of year in which a certain condition is met.

My current approach is the following (APHRO being the raster brick), but returns the error message below:

enter code here
r <- raster(ncol=40, nrow=20)
r[] <- rnorm(n=ncell(r))
APHRO <- brick(x=c(r, r*2, r))    
NewLayer <- calc(APHRO, fun=FindOnsetDate(APHRO))

Returning this error:

 Error in .local(x, ...) : not a valid subset 

And the function being parsed:

FindOnsetDate <- function (s) {
  x=0  
  repeat {
    x+1
    if(s[[x]] >= 20 | s[[x]] + s[[x+1]] >= 20 & ChkFalseOnset() == FALSE)
    {break}
  }
  return(x);
}

With the function for the 3rd condition being:

ChkFalseOnset <- function (x) {



  for (i in 0:13){

    if (sum(APHRO[[x+i:x+i+7]]) >= 5)
    {return(FALSE); break}
    return(TRUE)  
  }
}

Thank you in advance!!!!

And please let me know if I should provide more information - tried to keep it parsimonious.

AntonUrf
  • 21
  • 6

2 Answers2

0

I would suggest a basic two-step process. With a 365-days example:

set.seed(123)

r <- raster(ncol=40, nrow=20)

r_list <- list()

for(i in 1:365){
  r_list[[i]] <- setValues(r,rnorm(n=ncell(r),mean = 10,sd = 5))
}

APHRO <- brick(r_list)

Use a basic logic test for each iteration:

r_list2 <- list()

for(i in 1:365){
  if(i != 365){
    r_list2[[i]] <- APHRO[[i]] >= 20 | APHRO[[i]] + APHRO[[i+1]] >= 20
  }else{
    r_list2[[i]] <- APHRO[[i]] >= 20
  }
}

Compute sum by year:

NewLayer <- calc(brick(r_list2), fun=sum)

plot(NewLayer)

enter image description here

aldo_tapia
  • 1,153
  • 16
  • 27
0

The problem is that your function is no good:

FindOnsetDate <- function (s) {
   x=0  
   repeat {
     x+1
     if(s[[x]] >= 20 | s[[x]] + s[[x+1]] >= 20)
     {break}
   }
   return(x);
 }
 FindOnsetDate(1:100)
 #Error in s[[x]] : 
 # attempt to select less than one element in get1index <real>

Perhaps something like this:

FindOnsetDate <- function (s) {
   j <- s + c(s[-1], 0)
   sum(j > 20 | s > 20)
   # if all values are positive, just do sum(j > 20)
}
FindOnsetDate(1:20)
#10

This works now:

 r <- calc(APHRO, FindOnsetDate)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Dear RobertH, @aldo_tapia, thank you for your help. It would be great if you could elaborate a little on you solution. I can't make it work on my specific problem - I am new to R and wrapping my head around it. More context: The data set holds daily precipitation data for one year over the Asian region. I am trying to determine the monsoon onset defined as the following: First day of 2 consecutive days that receive more than 20mm rainfall and do not experience a dry period of 7 days < 5mm in the 20 days thereafter. I tried the latter in the ChkFalseOnset - see edit. Thank you in advance! – AntonUrf Dec 17 '17 at 06:43
  • Your problem is really about how to write a somewhat complicated function in R. I suggest posting that as a new question, where you leave out all the spatial data, and present a vector with rainfall data and a first attempt at a function to identify the onset of the monsoon. – Robert Hijmans Dec 17 '17 at 17:04
  • Thanks Robert, will do as suggested. – AntonUrf Dec 18 '17 at 06:56
  • Late reply - but thanks for your guidance. After revisiting fundamentals of R and Raster data types it worked wonderfully and I could adjust as needed. – AntonUrf Feb 20 '18 at 09:09