0

I am trying to simultaneously discretize 5 stacked rasters in R into quartile values. I've written the following loop to do so, but it doesn't seem to be working correctly. In the code, "stack.disc" is the stack and "quartiles" is a 5 column data frame with the 5 rasters as columns and their quartile values listed in the rows.

for(i in 1:ncol(quartiles)) {
  for(j in 1:length(stack.disc@layers[[i]])) {
    if(isTRUE(stack.disc@layers[[i]]@data@values[j] >= quartiles[1,i] &
                stack.disc@layers[[i]]@data@values[j] < quartiles[2,i])) {
        stack.disc@layers[[i]]@data@values[j] = 1 
    }
    if(isTRUE(stack.disc@layers[[i]]@data@values[j] >= quartiles[2,i] &
                stack.disc@layers[[i]]@data@values[j] < quartiles[3,i])) {
        stack.disc@layers[[i]]@data@values[j] = 2
    }
    if(isTRUE(stack.disc@layers[[i]]@data@values[j] >= quartiles[3,i] &
                stack.disc@layers[[i]]@data@values[j] < quartiles[4,i])) {
        stack.disc@layers[[i]]@data@values[j] = 3
    }
    if(isTRUE(stack.disc@layers[[i]]@data@values[j] >= quartiles[4,i] &
                stack.disc@layers[[i]]@data@values[j] <= quartiles[5,i])) {
        stack.disc@layers[[i]]@data@values[j] = 4
    }
  }
}

The code runs but only on the first 3 rasters in the stack. Any ideas why it isn't working for the last 2?

Thank you!

Jay Schoen
  • 3
  • 2
  • Maybe you should use one statement with many `else if ...` instead of five separate `if ...` statements – mob Mar 11 '20 at 00:02

1 Answers1

1

You can use cut to reclassify rasters (it has a raster* method).

s2 = stack(lapply(names(stack.disc), function(n) 
                  cut(stack.disc[[n]], breaks=quartiles[,n])))

Or, slightly shorter version using stackApply

stackApply(stack.disc, 1:5, function(r, na.rm) cut(r, breaks=quartiles[,names(r)]))

Some reproducible data to demonstrate:

stack.disc = stack(lapply(1:5, function(i) raster(matrix(rnorm(25, i),5,5))))
quartiles = t(quantile(stack.disc))
#         layer.1  layer.2  layer.3  layer.4  layer.5
# 0%   -1.0937068 0.498138 1.142862 2.229032 3.078026
# 25%   0.5171343 1.799564 2.496730 3.108751 4.484395
# 50%   1.1293477 2.099162 2.896269 3.818627 4.939167
# 75%   1.6348539 2.481976 3.615938 4.693733 5.314098
# 100%  2.4405652 3.511051 4.886970 5.456095 6.929452
dww
  • 30,425
  • 5
  • 68
  • 111
  • Much simpler and cleaner, thank you!. I found an error earlier in my code preventing two of the rasters from being processed properly (for both my code and yours), but it worked after catching that. Much appreciated. – Jay Schoen Mar 11 '20 at 14:36
  • You're welcome. For future reference, if you ever find yourself trying to access raster slots directly (the `@` operator) its a red flag that there's almost certainly a simpler cleaner approach available using other functions from the `raster` package. – dww Mar 11 '20 at 15:53