0

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

4

1 回答 1

1

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

如您所愿，我曾经st_apply()计算 class 的两个栅格的波段的所有成对组合的乘积stars。

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

  1. 输入：
  • 两个stars物体或两个stars_proxy物体
  • 两个栅格之间的波段数可以相同或不同
  1. 输出：
  • 一个stars具有名为“bandsProducts”维度的对象，其中包含两个栅格波段之间的所有成对产品组合。
  • 每个产品都有一个形式为 r2bX*r1bY 的名称（即一个值 - 参见下面的 Reprex），其中 X 和 Y 分别是栅格 2 和 1 的波段编号。

代表：

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 包于 2021-09-24 创建(v2.0.1)

于 2021-09-23T22:53:00.470 回答