indexing multi-layer rasters

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

indexing multi-layer rasters

Ben Tupper
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




--
Ben Tupper
Bigelow Laboratory for Ocean Science
West Boothbay Harbor, Maine
http://www.bigelow.org/
https://eco.bigelow.org

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: indexing multi-layer rasters

Bede-Fazekas Ákos

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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: indexing multi-layer rasters

Vijay Lulla
Hmmm... maybe the below might be easier to remember (at least it is for me)
then:

logo[[which(c(TRUE, TRUE, TRUE))]]
logo[[which(c(TRUE, FALSE, TRUE))]]

Thanks Akos!

On Fri, Jan 17, 2020 at 7:38 AM Bede-Fazekas Ákos <[hidden email]>
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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


--
Vijay Lulla

Assistant Professor,
Dept. of Geography, IUPUI
425 University Blvd, CA-207C.
Indianapolis, IN-46202
[hidden email]
ORCID: https://orcid.org/0000-0002-0823-2522
Online: http://vijaylulla.com | https://github.com/vlulla

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: indexing multi-layer rasters

Ben Tupper
In reply to this post by Bede-Fazekas Ákos
Ahhhh, the lightbulb just went off!  That makes perfect sense.

Thanks!
Ben

On Fri, Jan 17, 2020 at 10:45 AM Bede-Fazekas Ákos <[hidden email]> 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 <[hidden email]> 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
> >> [hidden email]
> >> 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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo