How can I create a multiband raster from singleband raster values?

Question

Counter-intuitive to all the other questions about splitting raster bands, I would like to split my single-band raster into a multiband raster stack based on the values given in each pixel (range from 1 to 11). Therefore, I would like to conditionally extract pixels based on their values and place them in an individual band. Basically:

band1 <- (if rastervalue == 1)
band2 <- (if rastervalue == 2)
.
.
.
band11 <- (if rastervalue == 11)

stacked <- stack(band1,band2,...band11)
writeRaster("path.tif")

Here is the overview of the data:

class      : RasterLayer 
dimensions : 25658, 15349, 393824642  (nrow, ncol, ncell)
resolution : 249.9999, 249.9999  (x, y)
extent     : -15191755, -11354506, 8399516, 14814014  (xmin, xmax, ymin, ymax)
crs        : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs 
values     : 1, 11  (min, max)

I have tried to do the following:

library(raster)
data <- raster("path.tif")
band1 = data[data != 1] = NA
band2 = data[data != 2] = NA
...
band11 = data[data != 11] = NA
stacked <- stack(band1,band2,...band11)

But, I get raster layers of NA for all the bands I create How can I do that in R?

Answer 1

Test raster with numbers 1 to 6 twice:

> r = raster(matrix(1:6, 3,4))

Loop over unique values in r, making a raster from each value, then feed to stack via do.call to make a stack:

> s = do.call(stack, lapply(unique(r[]), function(v){r==v}))
> plot(s)

I'm not sure if you want layers of TRUE/FALSE (1/0) or if you want to put the v values in the layer and NA elsewhere. You can make the function anything you like.

This does zeroes and v values:

s = do.call(stack, lapply(unique(r[]), function(v){(r==v)*v}))

This does NA and v values (but I can't help thinking there's an easier way):

s = do.call(stack, lapply(unique(r[]), function(v){t=r;t[]=NA;t[r==v]=v;t}))

Inspired by my comment on local scoped data, this would be the "easier" way - r is scoped within the function so doesn't change the global r here:

s = do.call(stack, lapply(unique(r[]), function(v){r[r!=v]=NA;r}))

If you know in advance the values then replace unique(r[]) with a vector of the values, otherwise R might take a while to find them. Your raster is quite large so this might all take a while anyway.