Products of all pairwise combinations of bands in two rasters with R stars package

Question

I have two multiband rasters of class stars. They have the same resolution and extent in their first two dimensions (x and y). Each raster has multiple bands. I would like to take all pairwise combinations of bands from each of the rasters and find the product of each of those combinations. Is there a way to do this with a function like outer() or possibly st_apply(), without having to use nested for-loops?

Answer 1

Hoping that it is not too late and that my answer will still be useful for you @qdread, I suggest the following solution (see Reprex below).

As you wished, I used st_apply() to compute the products of all pairwise combinations of bands of the two rasters of class stars.

For your convenience, I have built a function (i.e. named crossBandsProducts()) that wraps the whole process. This function has the following features:

1. Input:

• Two stars objects or two stars_proxy objects
• the number of bands can be the same or different between the two rasters

2. Output:

• One stars object with a dimension named "bandsProducts" containing all pairwise combinations of products between the bands of the two rasters.
• Each product has a name (i.e. a value - see Reprex below) of the form r2bX*r1bY where X and Y are respectively the band numbers of rasters 2 and 1.

REPREX:

library(stars)
#> Le chargement a nécessité le package : abind
#> Le chargement a nécessité le package : sf
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

# 1. Importing two stars objects with 6 and 3 bands respectively
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
tif1 <- read_stars(tif, proxy = FALSE)
tif2 <- read_stars(tif, proxy = FALSE, RasterIO = list(bands = c(1, 3, 4)))


# 2. Building the 'crossBandsProducts' function
crossBandsProducts <- function(r1, r2) {
  products <-
    st_apply(r2, 3, function(x)
      x * as.data.frame(split(r1, "band"))[, -grep("x|y", colnames(as.data.frame(split(r1, "band"))))])
  
  if (class(r1)[1] == "stars_proxy"){
    products <- st_as_stars(products)
  }
  
  products <- as.data.frame(split(products[[1]], "band"))
  
  colnames(products) <-
    paste0(rep(paste0("r2b", 1:dim(r2)["band"]), each = dim(r1)["band"]), 
           rep(paste0("*r1b", 1:dim(r1)["band"]), times = dim(r2)["band"]))
  
  products <-
    cbind(as.data.frame(split(r1, "band"))[, grep("x|y", colnames(as.data.frame(split(r1, "band"))))], products)
  
  products <- st_as_stars(products, dims = c("x", "y"))
  
  st_dimensions(products) <- st_dimensions(r1)[c("x", "y")]
  
  products <- st_set_dimensions(merge(products),
                                names = c("x", "y", "bandsProducts"))
  
  return(products)
}


# 3. Use of the function 'crossBandProducts'
(products <- crossBandsProducts(r1=tif1, r2=tif2))
#> stars object with 3 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>    Min. 1st Qu. Median     Mean 3rd Qu.  Max.
#> X  2209    4225   5776 6176.508    7569 65025
#> dimension(s):
#>               from  to  offset delta                       refsys point
#> x                1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE
#> y                1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE
#> bandsProducts    1  18      NA    NA                           NA    NA
#>                                values x/y
#> x                                NULL [x]
#> y                                NULL [y]
#> bandsProducts r2b1*r1b1,...,r2b3*r1b6
#>
#> NB: The dimension 'bandsProducts' has 18 values, which is consistent since the
#> rasters tif1 and tif2 have 6 and 3 bands respectively.


# 4. Example of extraction : two possibilities
# 4.1. via the number(s) of 'bandsProducts'
products[,,,4]
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>    Min. 1st Qu. Median     Mean 3rd Qu.  Max.
#> X   649    3904   4788 4528.267    5525 65025
#> dimension(s):
#>               from  to  offset delta                       refsys point
#> x                1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE
#> y                1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE
#> bandsProducts    4   4      NA    NA                           NA    NA
#>                  values x/y
#> x                  NULL [x]
#> y                  NULL [y]
#> bandsProducts r2b1*r1b4

# 4.2. via the name(s) (i.e. value(s)) of 'bandsProducts'
products[,,, values = c("r2b1*r1b4", "r2b3*r1b6")]
#> stars object with 3 dimensions and 1 attribute
#> attribute(s), summary of first 1e+05 cells:
#>    Min. 1st Qu. Median     Mean 3rd Qu.  Max.
#> X  2209    4225   5776 6176.508    7569 65025
#> dimension(s):
#>               from  to  offset delta                       refsys point
#> x                1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE
#> y                1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE
#> bandsProducts    1  18      NA    NA                           NA    NA
#>                                values x/y
#> x                                NULL [x]
#> y                                NULL [y]
#> bandsProducts r2b1*r1b1,...,r2b3*r1b6


# 5. Example of visualization
# 5.1. All pairwise combinations of bands products
plot(products, axes = TRUE, key.pos = NULL)
#> downsample set to c(2,2,1)

# 5.2. Selected pairwise combinations of bands products (selection by names/values)
plot(products[,,, values = c("r2b1*r1b4", "r2b3*r1b6")], axes = TRUE, key.pos = NULL)
#> downsample set to c(2,2,1)
#>
#> NB: Don't know why this second figure doesn't appear in the reprex, but in any
#> case it displays without any problem on my computer, so you shouldn't have any 
#> problem to display it when running this reprex locally.

^{Created on 2021-09-24 by the reprex package (v2.0.1)}