I have a large RasterStack (114 geotiff images) that I have successfully geoprocessed (masked, etc.) in R, but I am having difficulty getting it to apply a simple conditional statement to each raster (all raster layers are the same extent, resolution, and are co-registered). I want to set all pixel values that are less than 95% of each raster's max value to NA. For example, if a layer's max was 85, then pixel values < 80.75 = NA. Here's my code:
#Get max value from each raster layer
r_max <- maxValue(rstack)
#Set all values < 95% of max to NA
rstack[rstack < (r_max * 0.95)] = NA
When I run this code on the entire raster stack, I get "Error in value[j, ] : incorrect number of dimensions." However, if I run it on a smaller set (14 or so), it works exactly as it should. Because I have successfully executed a number of similar operations (other conditional statements, masking,etc.) on the entire stack without error, I am not sure why it's throwing this error now. Any ideas?
I apologize if this has been discussed before, but I was unable to find such a post. If it does exist, please point me in that direction.
Thanks
Here is a minimal reproducible example:
library(raster)
s <- stack(system.file("external/rlogo.grd", package="raster"))
cutoff <- maxValue(s) * .95
cutoff
#[1] 242.25 242.25 242.25
Now you can do:
s[s < cutoff] = NA
s
#class : RasterBrick
#dimensions : 77, 101, 7777, 3 (nrow, ncol, ncell, nlayers)
#resolution : 1, 1 (x, y)
#extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
#crs : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#source : memory
#names : red, green, blue
#min values : 243, 243, 243
#max values : 255, 255, 255
But there is a bug when the RasterStack is large (and needs to be written to file) --- and that is what you have stumbled upon. We can emulate that situation with rasterOptions(todisk=TRUE):
rasterOptions(todisk=TRUE)
s[s < cutoff] = NA
#Error in value[j, ] : incorrect number of dimensions
I will try to fix that. Here is workaround:
s <- stack(system.file("external/rlogo.grd", package="raster"))
cutoff <- maxValue(s) * .95
x <- sapply(1:nlayers(s), function(i) reclassify(s[[i]], cbind(-Inf, cutoff[i], NA)))
x <- stack(x)