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

【问题标题】：How to identify the polygons in which raster values were extracted with the stars R package?如何识别使用 stars R 包提取栅格值的多边形？
【发布时间】：2021-06-05 04:34:07
【问题描述】：

在上一个问题 (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-03-06 创建}

【问题讨论】：

标签： r raster r-stars

【解决方案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 年 9 月 21 日创建}

【讨论】：