1

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?

2
  • 1
    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!
    – Spacedman
    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!
    – Thomas
    Sep 23 at 13:31
2

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)

enter image description here

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}))

enter image description here

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.