如何使用带有 R 星包的多边形从栅格中提取值？

Question

使用stars 包，st_extract() 函数可以从定义位置的栅格中提取值。

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

tif <- system.file("tif/L7_ETMs.tif", package = "stars")
r <- read_stars(tif)
pnt <- st_sample(st_as_sfc(st_bbox(r)), 10)

st_extract(r[,,,1], pnt)
#> Simple feature collection with 10 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 288937.2 ymin: 9112173 xmax: 298589.9 ymax: 9120349
#> projected CRS:  UTM Zone 25, Southern Hemisphere
#>    L7_ETMs.tif                 geometry
#> 1           64 POINT (294613.4 9117565)
#> 2           72   POINT (295130 9117225)
#> 3           94 POINT (298589.9 9116806)
#> 4           86 POINT (296430.2 9112864)
#> 5           87 POINT (297481.9 9115176)
#> 6          110 POINT (288937.2 9112173)
#> 7           63 POINT (290966.6 9116890)
#> 8           84 POINT (295595.5 9116938)
#> 9           73 POINT (291047.1 9120349)
#> 10          65 POINT (294525.2 9117110)

我想做的是在这些点周围使用一个缓冲区并提取，比如说，缓冲区内的mean 值。 创建缓冲区

poly <- st_buffer(pnt, dist = 100)

现在我们有了多边形

poly
#> Geometry set for 10 features 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 288837.2 ymin: 9112073 xmax: 298689.9 ymax: 9120449
#> projected CRS:  UTM Zone 25, Southern Hemisphere
#> First 5 geometries:
#> POLYGON ((294713.4 9117565, 294713.3 9117560, 2...
#> POLYGON ((295230 9117225, 295229.8 9117220, 295...
#> POLYGON ((298689.9 9116806, 298689.8 9116800, 2...
#> POLYGON ((296530.2 9112864, 296530.1 9112859, 2...
#> POLYGON ((297581.9 9115176, 297581.8 9115171, 2...

问题就在这里，st_extract() 函数只使用点而不是多边形。

st_extract(r[,,,1], poly)
#> Error in st_extract.stars(r[, , , 1], poly): all(st_dimension(pts) == 0) is not TRUE

有没有办法提取多边形下的信息？

^{由reprex package (v1.0.0) 于 2021-02-19 创建}

Answer 1

这可以通过aggregate来完成：

library(stars)

tif = system.file("tif/L7_ETMs.tif", package = "stars")
r = read_stars(tif)
pnt = st_sample(st_as_sfc(st_bbox(r)), 10)
poly = st_buffer(pnt, dist = 100)

# Extract average value per polygon
x = aggregate(r, poly, mean)
st_as_sf(x)

## Simple feature collection with 10 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 289038 ymin: 9111186 xmax: 298491.2 ymax: 9120605
## projected CRS:  UTM Zone 25, Southern Hemisphere
##    L7_ETMs.tif.V1 L7_ETMs.tif.V2 L7_ETMs.tif.V3 L7_ETMs.tif.V4 L7_ETMs.tif.V5
## 1        87.80556       78.38889       87.44444       69.13889      124.05556
## 2        59.31579       43.94737       33.34211       70.76316       63.65789
## 3        78.33333       64.25641       62.56410       57.00000       70.79487
## 4        70.87179       57.89744       55.35897       63.94872       88.87179
## 5        89.51282       78.12821       86.00000       64.28205      111.48718
## 6        83.28205       67.46154       67.38462       51.38462       86.12821
## 7        80.27027       70.81081       72.59459       77.51351      103.78378
## 8        74.91892       60.75676       54.05405       85.86486       90.00000
## 9        68.56410       59.74359       55.10256       78.23077       94.41026
## 10       74.86486       60.91892       62.35135       58.91892      102.29730
##    L7_ETMs.tif.V6                       geometry
## 1        98.41667 POLYGON ((295003.7 9116093,...
## 2        31.55263 POLYGON ((290092.1 9119590,...
## 3        50.64103 POLYGON ((294767 9112633, 2...
## 4        61.38462 POLYGON ((289238 9114301, 2...
## 5        90.94872 POLYGON ((298491.2 9120505,...
## 6        69.41026 POLYGON ((289770 9111286, 2...
## 7        73.64865 POLYGON ((294775.7 9117676,...
## 8        57.78378 POLYGON ((294266.6 9113127,...
## 9        56.92308 POLYGON ((293838.6 9118091,...
## 10       77.51351 POLYGON ((290557.6 9114384,...

请记住，如果多边形之间存在重叠（与此示例不同），则每个栅格值仅在其落入的第一个多边形中“计数”一次（以符合aggregate 的普通行为）。