[R-sig-Geo] indexing multi-layer rasters
Ben Tupper
btupper @end|ng |rom b|ge|ow@org
Fri Jan 17 17:07:58 CET 2020
Ahhhh, the lightbulb just went off! That makes perfect sense.
Thanks!
Ben
On Fri, Jan 17, 2020 at 10:45 AM Bede-Fazekas Ákos <bfalevlist using gmail.com> wrote:
>
> 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 using 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
> >>>
> >>>
> >>>
> >>>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo using r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
>
--
Ben Tupper
Bigelow Laboratory for Ocean Science
West Boothbay Harbor, Maine
http://www.bigelow.org/
https://eco.bigelow.org
More information about the R-sig-Geo
mailing list