[GIS] Extract the cells of a raster based on a logical query R

extractrraster

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.

Best Answer

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

library(raster)

set.seed(123)

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_',
                                     'macrostachys_','orbignyana_',
                                     'stricta_','undata_'),'current')

currentStack_mask
## 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
mask(currentStack_mask[[1]],currentStack_mask[[1]],maskvalue=0)

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

Which is faster?

library(microbenchmark)

microbenchmark(first=mask(currentStack_mask[[1]],currentStack_mask[[1]],maskvalue=0),
               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){
  assign(names(currentStack_mask[[i]]),mask(currentStack_mask[[i]],currentStack_mask[[i]],maskvalue=0))
}
Related Question