indexing multi-layer rasters

Hello Ben,

I think you are absolutely right about "raster's implementation of `[[`
is different than base R". But I don't agree on your interpretation on
the logical indexing of rasters ('first logical index is used'). It is
not the first one. The logical vector is coerced to integer, and c(TRUE,
TRUE, TRUE) is treated as c(1, 1, 1), while c(TRUE, FALSE, TRUE) is
treated as c(1, 0, 1). This is why it cause a 'not a valid subset' error
(there is no sense of searching for the 0th layer of the RasterStack).
If the first logical index was used, c(TRUE, TRUE, TRUE) and c(TRUE,
FALSE, TRUE) would give exactly the same result, i.e a RasterLayer
created from the first layer of the RasterStack.

Have a nice weekend,
?kos


2020.01.17. 15:40 keltez?ssel, Ben Tupper ?rta:

Hi,

Thanks for this.  I think the root of my muddle is the mish-mash model
of a raster that I have in my head.  Depending upon what is most
convenient, I sometimes view a raster as a multi-dimensional array and
at other times as a list of single layer matrices.   If we step back
from logical indexing and use integers it is easier to identify
raster's slight variation on base R recursive indexing.  The example
below uses integer indices on a list and on a raster.

Back to logical indexing, in a way it makes perfect sense as just the
first logical index is used; just as if() does.  But what is different
is that it uses that first logical index for each element in the index
vector.  That's why logo[[c(TRUE, TRUE, TRUE)]] yields red.1, red.2
and red.3.

Thanks again and cheers,
Ben


library(raster)
x <- list(red = "R", green = "G", blue = "B")
logo <- raster::brick(system.file("external/rlogo.grd", package="raster"))

# `[` index a list, get a list back (with NULL for not found)
x[c(1,3)]
# $red
# [1] "R"
#
# $green
# [1] "G"

# `[` index a raster, get a matrix back (or vector for single layer)
logo[c(1,3)]
#      red green blue
# [1,] 255   255  255
# [2,] 255   255  255

# `[[` recursive index of a list fails
x[[c(1,3)]]
# Error in x[[c(1, 3)]] : subscript out of bounds

# `[[` index of a raster yields raster
# so raster's implementation of `[[` is different than base R
logo[[c(1,3)]]
# class      : RasterStack
# dimensions : 77, 101, 7777, 2  (nrow, ncol, ncell, nlayers)
# resolution : 1, 1  (x, y)
# extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
# crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# names      : red, blue
# min values :   0,    0
# max values : 255,  255

On Fri, Jan 17, 2020 at 7:38 AM Bede-Fazekas ?kos <bfalevlist at gmail.com> wrote:

Dear Ben,
Although I cannot answer your question on why logical subsetting was not
implemented in package raster, there is a very easy workaround:
logo[[(1:nlayers(logo))[c(TRUE, TRUE, TRUE)]]]
logo[[(1:nlayers(logo))[c(TRUE, FALSE, TRUE)]]]

Also note that in case of lists '[[' does recursive indexing, and this
type of logical indexing you are asking about works only with '['.
HTH,
?kos Bede-Fazekas
Hungarian Academy of Sciences

2020.01.16. 17:50 keltez?ssel, Ben Tupper ?rta:

Hi,

I usually avoid using logical indexes with multilayer rasters (stacks
and bricks), but I have never understood why indexing ala `[[` with
logicals isn't supported.  Below is an example that shows how it
doesn't work like the typical indexing with logicals.  I'm not asking
to have logical indices considered (although it would be handy), but
instead I really just want to understand why it is the way it is.  I
scanned over the introductory vignette (https://rspatial.org/raster)
and issues (https://github.com/rspatial/raster/issues) but found
nothing there.

Thanks and cheers,
Ben

### START
library(raster)

logo <- raster::brick(system.file("external/rlogo.grd", package="raster"))

# red-red-red
logo[[c(TRUE, TRUE, TRUE)]]
# class      : RasterStack
# dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
# resolution : 1, 1  (x, y)
# extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
# crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# names      : red.1, red.2, red.3
# min values :     0,     0,     0
# max values :   255,   255,   255

# red-red
logo[[c(TRUE, TRUE)]]
# class      : RasterStack
# dimensions : 77, 101, 7777, 2  (nrow, ncol, ncell, nlayers)
# resolution : 1, 1  (x, y)
# extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
# crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# names      : red.1, red.2
# min values :     0,     0
# max values :   255,   255

# red
logo[[c(TRUE)]]
# class      : RasterLayer
# band       : 1  (of  3  bands)
# dimensions : 77, 101, 7777  (nrow, ncol, ncell)
# resolution : 1, 1  (x, y)
# extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
# crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# source     : /usr/lib64/R/library/raster/external/rlogo.grd
# names      : red
# values     : 0, 255  (min, max)

logo[[c(TRUE, FALSE, TRUE)]]
#Error in .local(x, ...) : not a valid subset


#sessionInfo()
# R version 3.5.1 (2018-07-02)
# Platform: x86_64-redhat-linux-gnu (64-bit)
# Running under: CentOS Linux 7 (Core)
#
# Matrix products: default
# BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
#
# locale:
#   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
# [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
LC_PAPER=en_US.UTF-8       LC_NAME=C
# [9] LC_ADDRESS=C               LC_TELEPHONE=C
LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base
#
# other attached packages:
#   [1] raster_3.0-7 sp_1.3-2
#
# loaded via a namespace (and not attached):
#   [1] compiler_3.5.1   rgdal_1.4-8      tools_3.5.1      yaml_2.2.0
       Rcpp_1.0.3       codetools_0.2-15
# [7] grid_3.5.1       lattice_0.20-35

### END