Thanks to Friede, I've made a function to calculate the distance of many points along a line. Posting here in case this is useful for future reference. https://gist.github.com/wpetry/bb85a1ec3c408b2ab5dae17bd1e7771c

v_distance_along <- function(points, line, dist_unit = "km") {
  # packages
  require(sf)
  require(sfnetworks)
  require(dplyr)
  require(units)
  # check inputs
  if (!inherits(points, "sf") && !inherits(points, "sfc")) {
    stop("'points' must be an sf or sfc object containing POINT or MULTIPOINT geometries.")
  }
  if (!inherits(line, "sf") && !inherits(line, "sfc")) {
    stop("'line' must be an sf or sfc object containing a LINESTRING or MULTILINESTRING geometry.")
  }
  if (!all(st_geometry_type(points) %in% c("POINT", "MULTIPOINT"))) {
    stop("The second argument must be POINT or MULTIPOINT geometries.")
  }
  if (!st_geometry_type(line) %in% c("LINESTRING", "MULTILINESTRING")) {
    stop("The first argument must be a LINESTRING or MULTILINESTRING.")
  }
  if (is.na(sf::st_crs(points)) | is.na(sf::st_crs(line))) {
    stop("Both 'points' and 'line' must have a defined coordinate reference system (see ?st_crs).")
  }
  if (sf::st_is_longlat(points) | sf::st_is_longlat(line)) {
    stop("")
  }
  if (sf::st_crs(points) != sf::st_crs(line)) {
    stop("'points' and 'line' must have the same coordinate reference system (see ?st_crs).")
  }
  if (!units::ud_are_convertible(dist_unit, "meter")) {
    stop("'dist_unit' must be a valid unit of length.")
  }
  line <- sf::st_cast(line, "LINESTRING")  # ensure single LINESTRING
  path <- sfnetworks::as_sfnetwork(sf::st_as_sf(sf::st_sfc(line$geometry)),
                                   directed = FALSE)
  near <- sf::st_nearest_points(points, line)
  snap <- suppressWarnings(sf::st_sfc(lapply(near, function(l) sf::st_cast(l, "POINT")),
                                      crs = sf::st_crs(line)))
  pathx <- sfnetworks::st_network_blend(path, snap)  # add snapped points to network
  dist <- pathx |>
    sfnetworks::activate("edges") |>
    sf::st_as_sf() |>
    dplyr::mutate(length = sf::st_length(x)) |>
    sf::st_drop_geometry() |>
    dplyr::mutate(dist = round(cumsum(units::set_units(length - length[1], dist_unit,
                                                       mode = "standard")), 1),
                  from = dplyr::case_when(  # re-order vertices, moving line end to last position
                    from == 1L ~ 1L,
                    from == 2L ~ max(from),
                    from >= 3L ~ from - 1L,
                  )) |>
    dplyr::select(from, dist)
  return(dist)
}