2

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)

1 Answer 1

3

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).

3
  • Are there any alternatives to get this to run faster? Have a very high resolution DEM (2 meter horizontal res.) and want to extract averages over a large number of polygons.. Aggregate seems to run very very slow, although faster than using raster and raster::extract.
    – Heymans
    Commented Jul 3, 2021 at 22:21
  • 1
    The fastest option I'm aware of is the command line tool named exactextract (not the R package, which I found to be slower) github.com/isciences/exactextract Commented Jul 4, 2021 at 6:25
  • 3
    Thanks for the resource. I also found today that stars has a new, faster alternative by using st_extract but they are still an order of magnitude slower than exactextractr. See: github.com/r-spatial/stars/issues/421
    – Heymans
    Commented Jul 4, 2021 at 16:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.