带有 R 星包的两个栅格中所有成对波段组合的乘积

Question

我有两个 stars 类的多波段栅格。它们在前两个维度（x 和 y）中具有相同的分辨率和范围。每个栅格都有多个波段。我想从每个栅格中获取所有成对的波段组合，并找到每个组合的乘积。有没有办法用outer() 或st_apply() 之类的函数来做到这一点，而不必使用嵌套的for循环？

Answer 1

希望现在还为时不晚，并且我的回答仍然对您有用@qdread，我建议以下解决方案（请参阅下面的 Reprex）。

如您所愿，我使用st_apply() 计算了stars 类的两个栅格的波段的所有成对组合的乘积。

为了您的方便，我构建了一个函数（即命名为crossBandsProducts()）来包装整个过程。 该功能具有以下特点：

1. 输入：

• 两个stars 对象或两个stars_proxy 对象
• 两个栅格之间的波段数可以相同或不同

1. 输出：

• 一个 stars 对象，其维度名为“bandsProducts”，包含两个栅格波段之间的所有成对产品组合。
• 每个产品都有一个 r2bX*r1bY 形式的名称（即一个值 - 参见下面的 Reprex），其中 X 和 Y 分别是栅格 2 和 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.

^{由reprex package (v2.0.1) 于 2021 年 9 月 24 日创建}