如何识别使用stars R包提取栅格值的多边形？

如何解决如何识别使用stars R包提取栅格值的多边形？

跟进上一个问题 (on stackoverflow)，我试图了解使用多边形子集如何与 stars R 包配合使用。以下代码打开一个光栅文件并将其裁剪为较小的尺寸。

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.0,GDAL 3.2.1,PROJ 7.2.1
library(sf)
library(ggplot2)

tif <- system.file("tif/L7_ETMs.tif",package = "stars")
r <- read_stars(tif)[,1]

r <- r %>%
  st_crop(st_bBox(c(
    xmin = 294000,xmax = 294500,ymin = 9110800,ymax = 9111200
  ),crs = st_crs(r)))

现在我在这个网格上随机选取 4 个点。

set.seed(123)

pts <- st_sample(st_as_sfc(st_bBox(r)),4)

plot(r,key.pos = NULL,reset = FALSE)
plot(pts,add = TRUE,pch = 21,cex = 2,bg = "red",col = "red")

我将使用这四个点在每个点周围创建 30 米的缓冲区。

poly <- st_buffer(pts,dist = 30)

然后我可以按如下方式提取缓冲区下的值（创建一个 stars 对象）。

r[poly]
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>   L7_ETMs.tif   
#>  Min.   :71.00  
#>  1st Qu.:72.00  
#>  Median :74.50  
#>  Mean   :75.36  
#>  3rd Qu.:77.75  
#>  Max.   :85.00  
#>  NA's   :241    
#> dimension(s):
#>      from  to  offset delta                       refsys point values x/y
#> x     184 200  288776  28.5 UTM Zone 25,Southern Hem... FALSE   NULL [x]
#> y     336 350 9120761 -28.5 UTM Zone 25,Southern Hem... FALSE   NULL [y]
#> band    1   1      NA    NA                           NA    NA   NULL

使用 st_as_sf()，我可以将结果转换为多边形。

sf_poly <- st_as_sf(r[poly])
sf_poly
#> Simple feature collection with 14 features and 1 field
#> geometry type:  polyGON
#> dimension:      XY
#> bBox:           xmin: 293991.8 ymin: 9110786 xmax: 294476.3 ymax: 9111213
#> projected CRS:  UTM Zone 25,Southern Hemisphere
#> First 10 features:
#>    V1                       geometry
#> 1  80 polyGON ((294105.8 9111213,...
#> 2  85 polyGON ((294134.3 9111213,...
#> 3  79 polyGON ((294105.8 9111185,...
#> 4  71 polyGON ((294134.3 9111185,...
#> 5  78 polyGON ((294419.3 9111185,...
#> 6  73 polyGON ((294447.8 9111185,...
#> 7  77 polyGON ((294419.3 9111156,...
#> 8  72 polyGON ((294162.8 9111042,...
#> 9  72 polyGON ((294191.3 9111042,...
#> 10 76 polyGON ((294162.8 9111014,...

我们可以看到已经提取了 14 个像素。

ggplot() +
  geom_sf(data = sf_poly) +
  geom_sf(data = st_sfc(poly),fill = NA,color = "red") +
  theme_minimal()

我要问的问题是如何找出每个像素与哪个缓冲区相关联。例如，一个介于 1 和 4 之间的 id。

^{由 reprex package (v1.0.0) 于 2021 年 3 月 6 日创建}