跟进上一个问题(在 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 包于 2021-03-06 创建(v1.0.0)