I have a raster Brick which represents the distribution models of 7 palm species, named currentStack_mask which looks like this.

enter image description here

As you can see I have 7 species and each of their rasters are represented with 0 and 1 values.

Basically what I want to do is (for each species) to extract all the cells that have a value of 1 and create another raster with those cells and of course to do it in R because I want to keep a track of what I am doing and also is because is faster and don't have to deal with all the intermediate files.

The equivalent function in Arcgis of what I want to do is Spatial Analyst Tools -> Extraction -> Extract by Attributes, which basically extract the cells of a raster based on a logical query, which in this case is that the cell value is 1.

I have tried with extract() function of the Raster package but this function extract the values not the cells.

Can anybody help me?... I am sure there is a short way to do this.

  • Have you looked at the cellNumber argument for the extract function? You can also use getValues and use which on the output to get cell numbers that match your which criteria. I don't use raster bricks often enough to pull together an example.
    – JMT2080AD
    Oct 31, 2017 at 19:29

1 Answer 1


I'm thinking in two different ways to achieve this. First, I'll recreate your data:



r <- raster()

rlist <- list()

for(i in 1:7){
  rlist[[i]] <- setValues(r,sample(x=c(0,1),size=ncell(r),replace = T))

currentStack_mask <- stack(rlist)
names(currentStack_mask) <- paste0(c('cuneate_','deversa_','interrupta_',

## class       : RasterStack 
## dimensions  : 180, 360, 64800, 7  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## names       : cuneate_current, deversa_current, interrupta_current, macrostachys_current, orbignyana_current, stricta_current, undata_current 
## min values  :               0,               0,                  0,                    0,                  0,               0,              0 
## max values  :               1,               1,                  1,                    1,                  1,               1,              1 

I expose here two approachs:

# one approach

# second approach
currentStack_mask[[1]][currentStack_mask[[1]]==0] <- NA

Which is faster?


               second=currentStack_mask[[1]][currentStack_mask[[1]]==0] <- NA)
## Unit: milliseconds
##    expr      min        lq      mean   median       uq       max neval cld
##   first  4.59380  5.313997  5.690036  5.42912  5.65046  9.855744   100  a 
##  second 13.69026 14.078171 15.307921 14.65290 16.40512 21.504191   100   b

Let's use the first one... You can save each layer to a list or create new objects based on layer name (or other name). Also, If you want to save it, jus simply add writeRaster():

# all layer to a list (you can do a stack after)
outputs <- list()

for(i in 1:7){
  outputs[[i]] <- mask(currentStack_mask[[i]],currentStack_mask[[i]],maskvalue=0)

# each layer to a new object

for(i in 1:7){
  • 1
    Awesome solution +1, I've always used method #2 for this, will have to try #1 next time for efficiency.
    – GISHuman
    Oct 31, 2017 at 21:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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