0

I am trying to find the number of consecutive days that a certain criteria is met for a stack of rasters. In this case I am trying to figure out if there are 10 consecutive days without rain for each 30 day period of the year.

Because the extensive amount of data I am using I need to come up with an efficient way to do this. ie trying to avoid using a loop through all raster cells.

I found one function that might make things easier but I can't figure out how to efficiently apply it to all raster cells. Here is an example of how it might be done easily. Here is a way to very quickly calculate the number of consecutive TRUEs.

z <- c(TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE) rle(z) max(rle(z)$lengths)>=3

Now the issue is how (if possible) to apply that to all cells efficiently. stackApply? calc? Not sure.

Here is an example raster stack to help:

   set.seed(42)
    require(raster)
    r1 <- raster(nrows=10, ncols=10)
    r2=r3=r4=r5=r1
    r1[]= runif(ncell(r1))
    r2[]= runif(ncell(r1))
    r3[]= runif(ncell(r1))
    r4[]= runif(ncell(r1))
    r5[]= runif(ncell(r1))
    rs=stack(r1,r2,r3,r4,r5)<.25
mmann1123
  • 5,031
  • 7
  • 41
  • 49
  • Why are you using `raster` objects instead of simple matrices (or 3D `arrays` if necessary)? Beyond that, are you separating the year into 12 nonintersecting 30-day blocks, or do you want to run a rolling test over each 30-day sequence? – Carl Witthoft Jul 08 '14 at 19:43
  • They are read in as .tif files, so keeping them as raster stacks seems natural (just a multidimentional array) but willing to switch. Eventually I am going to have it rolling over each 30-day period. – mmann1123 Jul 08 '14 at 20:45
  • I keep it as a raster but loop over cells using extract. – mdsumner Jul 08 '14 at 21:43
  • use `calc(rs, fun)` where fun is your function (that works on a vector values, representing a single cell – Robert Hijmans Jun 23 '15 at 06:01

1 Answers1

0

Although I am sure there is a better way to do it. Here is one possibility:

threeconsecutive = function(rasterin){
     s = rasterToPoints( rasterin )
     X = s[,'x']
     Y = s[,'y']
     s = s[,3:dim(s)[2]]
     s  = apply(s,MARGIN=1, function(x){max( rle(c(x))$lengths[rle(c(x) )$values==T], na.rm=T )>=3} )
     return(rasterFromXYZ( data.frame(x=X,y=Y,true=s) ))
 }
out = threeconsecutive(rs)
plot(out)
mmann1123
  • 5,031
  • 7
  • 41
  • 49
  • Hi could you please help me with a related question. http://stackoverflow.com/questions/40467967/calculating-number-of-three-consecutive-values-above-a-threshold-in-raster-stac – rar Nov 07 '16 at 15:35