Using the stars package, it is possible to the st_extract() function to extract values from a raster at defined locations.
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)
What I would like to do is to use a buffer around these points and extract, let’s say, the mean values inside a buffer.
Create buffers
poly <- st_buffer(pnt, dist = 100)
Now we have polygons
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...
The problem is here, the st_extract() function uses only points and not polygons.
st_extract(r[,,,1], poly)
#> Error in st_extract.stars(r[, , , 1], poly): all(st_dimension(pts) == 0) is not TRUE
Is there a way to extract information under polygons?
Created on 2021-02-19 by the reprex package (v1.0.0)
This can be done with 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,...
Keep in mind that if there is overlap between polygons (unlike in this example) then each raster value is only "counted" once, in the first polygon it falls in (to comply with the ordinary behavior of aggregate).