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

Question

跟进上一个问题 ()，我试图了解使用多边形进行子集化如何与 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 日创建

Answer 1

可能有点晚了...但我向您推荐下面的解决方案，希望它仍然有用。

您只需在 reprex 的末尾添加以下代码行，它就可以正常工作。

library(dplyr)

# convert 'poly' from sfc to sf class object
poly <- st_as_sf(poly)

# create id's for the buffers (id's are equal to rows number) 
poly <- mutate(poly, id = row_number())

# intersects 'poly' with 'sf_poly' (i.e. identification of the pixels 
# intersected by each buffer)
(st_intersects(poly, sf_poly))
#> Sparse geometry binary predicate list of length 4, where the predicate
#> was `intersects'
#>  1: 1, 2, 3, 4
#>  2: 12, 13, 14
#>  3: 8, 9, 10, 11
#>  4: 5, 6, 7

# visualization
ggplot() +
  geom_sf(data = sf_poly) +
  geom_sf(data = st_geometry(poly), fill = NA, color = "red") +
  geom_sf_label(aes(label = poly$id, geometry = poly$x)) +
  theme_minimal()

^{由 reprex package (v2.0.1)}

于 2021-09-21 创建

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

How to identify the polygons in which raster values were extracted with the stars R package?

r

raster

r-stars