How to identify the polygons in which raster values were extracted with the stars R package?

Following-up a previous question (on stackoverflow), I am trying to understand how subsetting using polygons works with the stars R package. The following code opens a raster file and crops it to a smaller dimension.

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

tif <- system.file("tif/L7_ETMs.tif", package = "stars")
r <- read_stars(tif)[, , , 1]

r <- r %>%
  st_crop(st_bbox(c(
    xmin = 294000,
    xmax = 294500,
    ymin = 9110800,
    ymax = 9111200
  ), crs = st_crs(r)))

Now I randomly pick 4 points on this grid.

set.seed(123)

pts <- st_sample(st_as_sfc(st_bbox(r)), 4)

plot(r, key.pos = NULL, reset = FALSE)
plot(pts, add = TRUE, pch = 21, cex = 2, bg = "red", col = "red")

I will use these four points to create 30 meters buffers around each of them.

poly <- st_buffer(pts, dist = 30)

I can then extract the values under the buffers as follows (which create a stars object).

r[poly]
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>   L7_ETMs.tif   
#>  Min.   :71.00  
#>  1st Qu.:72.00  
#>  Median :74.50  
#>  Mean   :75.36  
#>  3rd Qu.:77.75  
#>  Max.   :85.00  
#>  NA's   :241    
#> dimension(s):
#>      from  to  offset delta                       refsys point values x/y
#> x     184 200  288776  28.5 UTM Zone 25, Southern Hem... FALSE   NULL [x]
#> y     336 350 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE   NULL [y]
#> band    1   1      NA    NA                           NA    NA   NULL

Using st_as_sf(), I can convert the results into polygons.

sf_poly <- st_as_sf(r[poly])
sf_poly
#> Simple feature collection with 14 features and 1 field
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 293991.8 ymin: 9110786 xmax: 294476.3 ymax: 9111213
#> projected CRS:  UTM Zone 25, Southern Hemisphere
#> First 10 features:
#>    V1                       geometry
#> 1  80 POLYGON ((294105.8 9111213,...
#> 2  85 POLYGON ((294134.3 9111213,...
#> 3  79 POLYGON ((294105.8 9111185,...
#> 4  71 POLYGON ((294134.3 9111185,...
#> 5  78 POLYGON ((294419.3 9111185,...
#> 6  73 POLYGON ((294447.8 9111185,...
#> 7  77 POLYGON ((294419.3 9111156,...
#> 8  72 POLYGON ((294162.8 9111042,...
#> 9  72 POLYGON ((294191.3 9111042,...
#> 10 76 POLYGON ((294162.8 9111014,...

We can see that there are 14 pixels that have been extracted.

ggplot() +
  geom_sf(data = sf_poly) +
  geom_sf(data = st_sfc(poly), fill = NA, color = "red") +
  theme_minimal()

The question I am asking is how can I find out to which buffer is associated each of these pixels. For example, an id between 1 and 4.

^{Created on 2021-03-06 by the reprex package (v1.0.0)}

Maybe a bit late... but I propose you the solution below and hope it will still be useful.

You just need to add the following lines of code at the end of your reprex and it should work.

library(dplyr)

# convert 'poly' from sfc to sf class object
poly <- st_as_sf(poly)

# create id's for the buffers (id's are equal to rows number) 
poly <- mutate(poly, id = row_number())

# intersects 'poly' with 'sf_poly' (i.e. identification of the pixels 
# intersected by each buffer)
(st_intersects(poly, sf_poly))
#> Sparse geometry binary predicate list of length 4, where the predicate
#> was `intersects'
#>  1: 1, 2, 3, 4
#>  2: 12, 13, 14
#>  3: 8, 9, 10, 11
#>  4: 5, 6, 7

# visualization
ggplot() +
  geom_sf(data = sf_poly) +
  geom_sf(data = st_geometry(poly), fill = NA, color = "red") +
  geom_sf_label(aes(label = poly$id, geometry = poly$x)) +
  theme_minimal()

^{Created on 2021-09-21 by the reprex package (v2.0.1)}