How to loop through polygons that overlap and randomly sample points from each in r

Is there a way to randomly sample 5 points from individual polygons using st_sample if those polygons overlap? I want the sampled points to inherit the attributes from the polygons they are sampled from.

# create an empty list to store results
empty_list <- list()

# loop through each unique polygon and sample 5 points from each
for(i in used_buff_1k$unique_id[i]) {
  empty_list <- st_sample(used_buff_1k, size = rep(5, nrow(used_buff_1k)), type = "random")
}

This created a list of 50 points (5 from 10 polygons), but it doesnt inherit the attributes, it just produced a list of points. It looks like this may have produced what I wanted, but since it didnt loop through each unique polygon and sample individually, there seems to be no way for me to link up which random points go with each polygon.

If tidyverse methods are OK, we could use purrr::map() in mutate() to generate a list of 5 sample points for each polygon. When resulting list-column is unnested, we will have a sf object with polygons and points, rows for each sampled point and all attributes of the original polygon sf. We can now switch the geometry to points column and drop columns we don't need.

library(sf)
library(dplyr)
library(ggplot2)
library(tidyr)

# generate sample data, 5 overlapping circles + id attribute
set.seed(42)
c_df <- list(matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)) %>% 
  st_polygon() %>% 
  st_sample(size = 5) %>% 
  st_buffer(4) %>% 
  st_as_sf() %>% 
  st_set_geometry("geometry") %>% 
  mutate(id = row_number()) %>% 
  as_tibble() %>% st_as_sf()

c_df
#> Simple feature collection with 5 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -1.138605 ymin: -2.653334 xmax: 13.37075 ymax: 11.36588
#> CRS:           NA
#> # A tibble: 5 × 2
#>                                                                   geometry    id
#>                                                                  <POLYGON> <int>
#> 1 ((13.14806 5.190959, 13.14258 4.981616, 13.12615 4.772846, 13.09881 4.5…     1
#> 2 ((13.37075 7.365883, 13.36527 7.156539, 13.34884 6.947769, 13.32151 6.7…     2
#> 3 ((6.861395 1.346666, 6.855913 1.137322, 6.839483 0.9285521, 6.812149 0.…     3
#> 4 ((12.30448 6.569923, 12.29899 6.360579, 12.28256 6.151809, 12.25523 5.9…     4
#> 5 ((10.41746 7.050648, 10.41197 6.841304, 10.39554 6.632534, 10.36821 6.4…     5

# sample 5 points from each polygon , 5 × 3 tibble, points are stored in a list-column 
# unnest sample_points list-column, resulting tibble is 25 × 3
# switch geometry column; resulting sf object includes poly_id with polygon id's
sampled_points <- c_df %>% 
  mutate(sample_points = purrr::map(geometry, st_sample, 5)) %>% 
  unnest_longer(sample_points) %>% 
  st_set_geometry("sample_points") %>% 
  select(poly_id = id)

Result:

sampled_points
#> Simple feature collection with 25 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 3.755624 ymin: 0.01408366 xmax: 13.28189 ymax: 9.854324
#> CRS:           NA
#> # A tibble: 25 × 2
#>    poly_id       sample_points
#>      <int>             <POINT>
#>  1       1 (8.809995 9.016771)
#>  2       1 (10.90096 2.130858)
#>  3       1 (12.62544 4.990936)
#>  4       1 (7.191491 5.673621)
#>  5       1 (8.846403 8.423211)
#>  6       2  (13.28189 6.94164)
#>  7       2 (6.030255 9.266648)
#>  8       2 (9.484448 9.854324)
#>  9       2 (8.492382 6.470749)
#> 10       2 (12.61666 8.847241)
#> # ℹ 15 more rows

mutate(sampled_points, poly_id = as.factor(poly_id)) %>% 
  ggplot() +
  geom_sf(data = mutate(c_df, id = as.factor(id)), aes(fill = id), alpha = .5) +
  geom_sf_text(aes(label = poly_id), nudge_x = -.25, nudge_y = -.25) +
  geom_sf(aes(fill = poly_id), shape = 21, size = 3) +
  scale_fill_viridis_d() +
  scale_color_viridis_d() +
  theme_void()

<!-- -->

^{Created on 2023-06-04 with reprex v2.0.2}