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?
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:
stars
objects or two stars_proxy
objectsstars
object with a dimension named "bandsProducts" containing all pairwise combinations of products between the bands of the two rasters.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)