Extraer las celdas de un raster a partir de una consulta lógica R

Question

Tengo un raster Brick que representa los modelos de distribución de 7 especies de palmeras, denominadas currentStack_mask que tiene este aspecto.

Como puedes ver tengo 7 especies y cada uno de sus raster están representados con valores 0 y 1.

Básicamente lo que quiero hacer es (para cada especie) extraer todas las celdas que tienen un valor de 1 y crear otro raster con esas celdas y por supuesto hacerlo en R porque quiero llevar un control de lo que estoy haciendo y además es porque es más rápido y no tengo que lidiar con todos los archivos intermedios.

La función equivalente en Arcgis de lo que quiero hacer es Spatial Analyst Tools -> Extraction -> Extract by Attributes, que básicamente extrae las celdas de un raster basándose en una consulta lógica, que en este caso es que el valor de la celda es 1.

He probado con extract() del paquete Raster, pero esta función extrae los valores, no las celdas.

¿Alguien puede ayudarme? Estoy seguro de que hay una manera corta de hacer esto.

Answer 1

Estoy pensando en dos formas diferentes de conseguirlo. En primer lugar, voy a recrear sus datos:

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 

Expongo aquí dos enfoques:

# one approach
mask(currentStack_mask[[1]],currentStack_mask[[1]],maskvalue=0)

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

¿Cuál es más rápido?

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

Usemos la primera... Puede guardar cada capa en una lista o crear nuevos objetos basados en el nombre de la capa (u otro nombre). Además, si desea guardarlo, sólo tiene que añadir 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))
}