Limitations with RasterStack math and conditional statements

Question

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

Robert Hijmans Robert Hijmans · Accepted Answer · 2020-03-26T04:06:41

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)

Limitations with RasterStack math and conditional statements

1 Answers