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)

• Charlodde: Thx a lot @Jeffrey Evans. This was exactly I was looking for. And everything is nicely explained that I could follow easily!