При профилировании вашей функции большую часть времени использует 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

^{Created on 2024-03-23 with 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
}

Вы можете использовать terra::perim, чтобы получить периметр многоугольников или длину линий.

library(terra)
v <- vect(system.file("ex/lux.shp", package = "terra"))
perim(v)
# [1] 117100.12  93477.28  84502.45  44919.14  85032.61  74708.05  57991.42  81203.93  74443.82  95564.74  80618.76 70810.67

Это работает как для плоских, так и для угловых (долгота/широта) координат.

Чтобы получить длину сегментов линий/многоугольников, вам необходимо сначала дезагрегировать их. Здесь показан первый многоугольник в v.

x <- disagg(v[1,], segments=TRUE)
p <- perim(x)
head(p)
#[1] 1383.51779  350.05402  580.19502   93.65138  333.32716  798.43439

sum(p)
#[1] 117100.1

И с вашими примерными данными вы получаете нужные цифры.

vect(lines) |> disagg(TRUE) |> perim() |> cumsum() |> head()
# [1] 1716.196 3290.379 4824.087 6745.759 7446.660 8303.283

Спасибо всем за ваши быстрые ответы; Мне нравятся все из них. В каждом предыдущем ответе был пример огромного увеличения скорости по сравнению с исходным вопросом. Я немного поигрался с кодом и нашел, насколько мне известно, наиболее эффективный подход. Это может быть полезно для дальнейшего использования. Вот тест с функцией Rust, написанный с помощью пакета расширения. Он превосходит {terra} примерно на 10%. 18 раз и основание R примерно на . 9 раз. Полученные векторы идентичны (см. check = TRUE).

  ## RUST
  # Wrap code from Josep Pueyo answer into RUST function
  code <- r"(
    use extendr_api::prelude::*;

    #[extendr]
    fn cartesian_distance(x1: f64, y1: f64, x2: f64, y2: f64) -> f64 {
      // Find max and min points
      let xmax = x1.max(x2);
      let xmin = x1.min(x2);
      let ymax = y1.max(y2);
      let ymin = y1.min(y2);

      // Calculate the distance between two points in 2D space
      let dx = xmax - xmin;
      let dy = ymax - ymin;
      (dx*dx + dy*dy).sqrt()
    }

    #[extendr]
    fn cartesian_distance_vector(x_vec: Vec<f64>, y_vec: Vec<f64>) -> Vec<f64> {

    // Pre-allocate result vector
    let mut cumulative_distances: Vec<f64> = Vec::with_capacity(x_vec.len());
    let n = x_vec.len() - 1;
    let mut cumulative_distance = 0.0;

    // Iterate over X-Y pairs
    for i in 0..n {

        let distance = cartesian_distance(x_vec[i], y_vec[i], x_vec[i+1], y_vec[i+1]);
        let distance = distance;

        // Accumulate cumulative distance
        cumulative_distance += distance;

        // Push cumulative distance to the result vector
        cumulative_distances.push(cumulative_distance);
    }

    cumulative_distances
}
)"
rextendr::rust_source(code = code)
#> ℹ build directory: 'C:/Users/TsyplenkovA/AppData/Local/Temp/RtmpquAnXb/file3d85de14e0f'
#> ✔ Writing 'C:/Users/TsyplenkovA/AppData/Local/Temp/RtmpquAnXb/file3d85de14e0f/target/extendr_wrappers.R'

## Base R
base_r <- function(X, Y){
  
  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)
  }
  
  n <- length(X) - 1
  
  # Pre-allocate a list 
  dists <- numeric(n)
  
  # Construct linestrings
  for (i in 1:n) {
    dists[[i]] <- cartesian_dist(X[i], Y[i], X[i+1], Y[i+1])
  }
  
  # Measure cumulative segment length
  cumsum(dists)
  
}

## Terra
library(terra)
#> terra 1.7.71
library(sf) # sf package is only needed for reading seine
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
lv <- vect(spData::seine[2, ])


## Benchmark
bench::mark(
  rust = cartesian_distance_vector(
    geom(lv)[, 3],
    geom(lv)[, 4]
  ),
  base_r = base_r(
    geom(lv)[, 3],
    geom(lv)[, 4]
  ),
  terra = {
    lv |>
      disagg(TRUE) |>
      perim() |>
      cumsum()
  },
  check = T,
  relative = T
)
#> # A tibble: 3 × 6
#>   expression   min median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <dbl>  <dbl>     <dbl>     <dbl>    <dbl>
#> 1 rust        1      1        16.9       10.1     4.71
#> 2 base_r      9.32   8.59      2.15      49.7     5.11
#> 3 terra      20.0   18.5       1          1       1

^{Created on 2024-03-26 with reprex v2.1.0}

Версия Rcpp даже быстрее, чем версия Rust:

library(Rcpp)
cppFunction('NumericVector cartesian_distance_cpp(NumericVector x_vec, NumericVector y_vec) {
  int n = x_vec.size() - 1;
  NumericVector cumulative_distances(n);
  double cumulative_distance = 0.0;

  for (int i = 0; i < n; i++) {
    double distance = sqrt(pow(x_vec[i] - x_vec[i+1], 2) + pow(y_vec[i] - y_vec[i+1], 2));
    cumulative_distance += distance;
    cumulative_distances[i] = cumulative_distance;
  }
  return cumulative_distances;
}')

На моей машине Rust в 1,6 раза медленнее, чем Rcpp... но учтите, что {terra} выигрывает с точки зрения памяти, поскольку требует в 7 раз меньше памяти, чем Rust, Rcpp или base.