st_explode

st_explode

library(sf)
library(mapview)
library(data.table)

lns = st_cast(trails, "LINESTRING")
# lns = lns[rep(seq_len(nrow(lns)), 1000),]
pts = st_cast(st_line_sample(lns, n = 1), "POINT")

st_explode = function(x) {
  if (!inherits(st_geometry(x), "sfc_LINESTRING")) 
    x = st_cast(x, "LINESTRING")
  crds = st_coordinates(x)[, 1:2]
  mat = lapply(seq(nrow(crds) - 1), function(i) {
    st_linestring(rbind(crds[i, ], crds[i + 1, ]))
  })
  st_sfc(mat, crs = st_crs(x))
}

extract_segment = function(ln, pt, radius) {
  ln_ex = st_explode(ln)
  ind = st_intersects(st_buffer(pt, radius), ln_ex)[[1]]
  st_sfc(ln_ex[[ind]], crs = st_crs(ln))
}

system.time({
  tst = st_as_sf(rbindlist(lapply(seq(nrow(lns)), function(i) {
    print(i)
    geom = extract_segment(lns[i, ], pts[i], 0.0001)
    as.data.table(st_sf(id = i, geom))
  })))
})