R: raster stack (time series) loop pixelwise and write layer number of value occurence into a single output layer

Question

I have a raster stack of gridded RFE 2.0 rainfall data (time series of one year = 365 layers). Now I wanted to identifiy on which date during the year the first rainfall occurs for each Pixel. The date is corresponding to the layer number (nlayers) = DOY (Day of the Year). The result should be a single output raster layer where the pixels are showing the layer number (=DOY) where the first rainfall occured. I am sure a loop would probably solve me problem, but I am really struggling. The occurence of rainfall is true if the value within a Pixel becomes > 0.

    library(raster)
    library(rgdal)

setwd("D:/*/2001") # Location of raster for e.g. year 2001
lst <- list.files(path=getwd(), pattern='.tif')
RFE2001 = stack(lst) 
RFE <- RFE2001

## Crop raster stack to extent of study area
e <- extent(10,21,6,15) # extent(xmin, xmax, ymin, ymax)
r <- crop(RFE, e)
crs(r) <- "+init=epsg:4326"

for (i in 1:nrow(r)) {
  for(j in 1:ncol(r)) {
    Value<-r[i,j]
    if (Value > 0) {
      s[i,j]<-Value
    } else {
   .....

Answer 1

You can pass a simple function, using "which" to return the desired julian day index, to the raster "calc" or "overlay" function. This will return a single raster layer with the first julian day of rain.

He we create some example data that approximates your problem with 20 raster layers in the stack.

library(raster)
r <- raster( xmn=10, xmx=21, ymn=6, ymx=15, resolution = 0.07, 
             crs = "+init=epsg:4326" )
    r[] <- runif(ncell(r), 0, 20)            
r <- stack(r)
  for(i in 1:20) {
    rx <- r[[1]]
    rx[] <- runif(ncell(rx), 0, 20)  
    r <- addLayer(r, rx)
  }  

Now we write a function "first.rain" that returns the index of the first raster with a value greater than "t" and pass it to "calc".

first.rain <- function(x, t = 1) { which(x > t)[1] }  
first <- calc(r, fun=first.rain)

To understand this a bit more we can create a simple example by creating a random vector.

( x <- runif(20, 0, 10) )

We apply "which" with the condition > 5. The function returns the position(s) of values that meet our condition.

which( x > 5 ) 

Note that it returns more than one value so, we can just return the first value using a bracket index.

which( x > 5 )[1] 

Now, what if we wanted the number of days it rained? We need to rethink our code because for a each pixel the number of value may be different. We can use something like "length" to solve this.

rain.days <- function(x, t = 1) { 
    return( length(x[x > t]) ) 
  }

rain.n <- calc(r, fun=rain.days)