0

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 {
   .....
4

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)
  • Thx a lot @Jeffrey Evans. This was exactly I was looking for. And everything is nicely explained that I could follow easily! – Charlodde Jan 26 '16 at 22:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.