Extract cell values at point locations

st_extract(x, ...)

# S3 method for stars
st_extract(
  x,
  at,
  ...,
  bilinear = FALSE,
  time_column = attr(at, "time_column") %||% attr(at, "time_col"),
  interpolate_time = bilinear,
  FUN = mean
)

Arguments

x

object of class stars or stars_proxy

...

passed on to aggregate.stars when geometries are not exclusively POINT geometries

at

object of class sf or sfc with geometries, or two-column matrix with points in rows, indicating where to extract x

bilinear

logical; use bilinear interpolation rather than nearest neighbour?

time_column

character or integer; name or index of a column with time or date values that will be matched to values of the dimension "time" in x, after which this dimension is reduced. This is useful to extract data cube values along a trajectory; see https://github.com/r-spatial/stars/issues/352 .

interpolate_time

logical; should time be interpolated? if FALSE, time instances are matched using the coinciding or the last preceding time in the data cube.

FUN

function used to aggregate pixel values when geometries of at intersect with more than one pixel

Value

if at is of class matrix, a matrix with extracted values is returned; otherwise: if x has more dimensions than only x and y (raster), an object of class stars with POINT geometries replacing x and y raster dimensions, if this is not the case, an object of sf with extracted values.

Details

points outside the raster are returned as NA values. For large sets of points for which extraction is needed, passing a matrix as to at may be much faster than passing an sf or sfc object.

Examples

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, pnt)
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>              Min. 1st Qu. Median     Mean 3rd Qu. Max.
#> L7_ETMs.tif    11      60     70 69.73333    86.5  121
#> dimension(s):
#>          from to offset delta                     refsys point
#> geometry    1 10     NA    NA SIRGAS 2000 / UTM zone 25S  TRUE
#> band        1  6     NA    NA                         NA    NA
#>                                                         values
#> geometry POINT (288862.6 9120298),...,POINT (298691.8 9111668)
#> band                                                      NULL
st_extract(r, pnt) %>% st_as_sf()
#> Simple feature collection with 10 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 288862.6 ymin: 9111668 xmax: 298691.8 ymax: 9120298
#> Projected CRS: SIRGAS 2000 / UTM zone 25S
#>    L7_ETMs.tif.V1 L7_ETMs.tif.V2 L7_ETMs.tif.V3 L7_ETMs.tif.V4 L7_ETMs.tif.V5
#> 1              72             64             68             71            107
#> 2              67             57             53             76             85
#> 3              71             59             59             79            108
#> 4              88             75             86             60            119
#> 5              98             93             76             14             13
#> 6              80             69             79             62            121
#> 7              88             74             81             60            119
#> 8              61             46             37             74             61
#> 9              75             66             69             64             98
#> 10            103             95             68             15             13
#>    L7_ETMs.tif.V6                 geometry
#> 1              66 POINT (288862.6 9120298)
#> 2              48   POINT (290373 9119475)
#> 3              60 POINT (293694.9 9116176)
#> 4              99 POINT (292121.2 9114392)
#> 5              12   POINT (297288 9112647)
#> 6              95 POINT (293561.9 9114499)
#> 7              98 POINT (293374.6 9112064)
#> 8              30 POINT (289304.6 9116860)
#> 9              69 POINT (296749.5 9120243)
#> 10             11 POINT (298691.8 9111668)
st_extract(r[,,,1], pnt)
#> Simple feature collection with 10 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 288862.6 ymin: 9111668 xmax: 298691.8 ymax: 9120298
#> Projected CRS: SIRGAS 2000 / UTM zone 25S
#>    L7_ETMs.tif                 geometry
#> 1           72 POINT (288862.6 9120298)
#> 2           67   POINT (290373 9119475)
#> 3           71 POINT (293694.9 9116176)
#> 4           88 POINT (292121.2 9114392)
#> 5           98   POINT (297288 9112647)
#> 6           80 POINT (293561.9 9114499)
#> 7           88 POINT (293374.6 9112064)
#> 8           61 POINT (289304.6 9116860)
#> 9           75 POINT (296749.5 9120243)
#> 10         103 POINT (298691.8 9111668)
st_extract(r, st_coordinates(pnt)) # "at" is a matrix: return a matrix
#>       [,1] [,2] [,3] [,4] [,5] [,6]
#>  [1,]   72   64   68   71  107   66
#>  [2,]   67   57   53   76   85   48
#>  [3,]   71   59   59   79  108   60
#>  [4,]   88   75   86   60  119   99
#>  [5,]   98   93   76   14   13   12
#>  [6,]   80   69   79   62  121   95
#>  [7,]   88   74   81   60  119   98
#>  [8,]   61   46   37   74   61   30
#>  [9,]   75   66   69   64   98   69
#> [10,]  103   95   68   15   13   11