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

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?

• Your code stomps on `data` when you do `data[data!=1]=NA` and then returns the value of that expression - not the new value of `data` - to store into `band1`. And the value of that expression is `NA`. You could do it by running in a `local` scope: `band1 = local({data[data!=1]=NA;data})` which leaves the main `data` object unchanged. Or see my answer! Sep 23 at 12:49
• Thanks. That is a good observation, I changed it and that also help. However, the answer from @Spacedman definitely did the job for me! Sep 23 at 13:31

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.

• Amazing! You saved me a lot of time - brilliant. Thank you. Sep 23 at 13:31