问题描述
使用 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[,poly)
#> Error in st_extract.stars(r[,poly): all(st_dimension(pts) == 0) is not TRUE
有没有办法提取多边形下的信息？
由 reprex package (v1.0.0) 于 2021 年 2 月 19 日创建
解决方法
这可以通过 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 ).