分析你的功能，大部分时间是geos_make_linestring.因此，假设您可以将数据转换为平面坐标(UTM)，则可以直接计算两点之间的距离来跳过该操作:

library(sf)
library(geos)
library(spData)
# Example
input <- seine[2, ] |> 
  st_transform(23032)

# My foo
cumulative_length <- function(input) {
  
  # Save CRS
  crs <- sf::st_crs(input)
  
  # Retrive coordinates
  lines_coo <- 
    sf::st_coordinates(input)
  
  # Count number of segments of linestring
  n <- nrow(lines_coo) - 1
  
  # Pre-allocate a list 
  lines_geos <- 
    vector(mode = "list", length = n)
  
  # Construct linestrings
  for (i in 1:n) {
    lines_geos[[i]] <- 
      geos::geos_make_linestring(lines_coo[i:(i+1),1], 
                                 lines_coo[i:(i+1),2], 
                                 crs = crs)
  }
  
  # Measure cumulative segment length
  lines_order <- 
    sapply(lines_geos, geos::geos_length) |> 
    append(0, 0) |> 
    cumsum()
  
  return(lines_order)
}

cartesian_dist <- function(x1, y1, x2, y2) {
  # √[(x2 − x1)2 + (y2 − y1)2]
  xmax <- max(x1, x2)
  xmin <- min(x1, x2)
  ymax <- max(y1, y2)
  ymin <- min(y1, y2)
  sqrt((xmax - xmin)^2 + (ymax - ymin)^2)
}

cumulative_length2 <- function(input) {
    
    # Retrive coordinates
    lines_coo <- sf::st_coordinates(input)
    
    # Count number of segments of linestring
    n <- nrow(lines_coo) - 1
    
    # Pre-allocate a list 
    dists <- numeric(n)
    
    # Construct linestrings
    for (i in 1:n) {
      dists[[i]] <- cartesian_dist(lines_coo[i, 1], lines_coo[i, 2], lines_coo[i+1, 1], lines_coo[i+1, 2])
    }
    
    # Measure cumulative segment length
    c(0, cumsum(dists))
    
}

microbenchmark::microbenchmark(
  cumulative_length(input),
  cumulative_length2(input),
  times = 100
)
#> Unit: milliseconds
#>                       expr     min      lq      mean   median       uq     max
#>   cumulative_length(input) 22.6733 23.5284 27.012571 24.42025 28.93890 53.6292
#>  cumulative_length2(input)  2.5730  2.7048  3.484514  2.82825  2.98925 21.8944
#>  neval
#>    100
#>    100

^{创建于2024—03—23，reprex v2.1.0}

如果无法转换为平面坐标，也可以使用WGS84投影并使用哈弗正弦公式计算点之间的距离，而不是笛卡尔距离:

haversine_dist <- function(lon1, lat1, lon2, lat2){
  R <- 6371e3 # metres
  phi1 <- lat1 * pi/180 # radians
  phi2 = lat2 * pi/180 # radians
  delta_phi = (lat2-lat1) * pi/180 # radians
  delta_lambda = (lon2-lon1) * pi/180
  
  a <- sin(delta_phi/2) * sin(delta_phi/2) + cos(phi1) * cos(phi2) * sin(delta_lambda/2) * sin(delta_lambda/2)
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  
  R * c
}