如何用sf/terra快速计算LINESTRING/MULTILINESTRING的累积长度？

需要计算LINESTRING或MULTILINESTRING从起点到各节点的累积距离，当前方法是拆分线段逐一测量，可行但耗时，求内置/更高效的实现方式。以下是基于spData塞纳河数据集的可复现示例，同时提供terra或rsgeo的解决方案。

原实现代码 #

library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library(geos)
library(spData)
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://nowosad.github.io/drat/', type='source')`

# Example
lines <- seine[2, ]

# 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)
    
  }

bench::mark(cumulative_length(lines))
#> # A tibble: 1 × 6
#>   expression                    min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>               <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 cumulative_length(lines)     13ms   13.7ms      72.9     244KB     109.

期望输出：

cumulative_length(lines) |> 
  head()
#> [1]    0.000 1716.196 3290.379 4824.087 6745.759 7446.660

高效解决方案 #

1. sf + geos 优化版 #

避免循环创建线段对象，直接计算相邻点的距离再累加，性能提升显著：

library(sf)
library(geos)
library(spData)

lines <- seine[2, ]

cumulative_length_fast <- function(input) {
  # 获取坐标矩阵
  coords <- st_coordinates(input)
  # 计算相邻点对的距离（自动适配CRS）
  seg_lengths <- geos::geos_distance(
    geos::geos_make_point(coords[-nrow(coords), 1], coords[-nrow(coords), 2], crs = st_crs(input)),
    geos::geos_make_point(coords[-1, 1], coords[-1, 2], crs = st_crs(input))
  )
  # 生成累积距离向量，开头补0
  c(0, cumsum(seg_lengths))
}

# 验证结果
cumulative_length_fast(lines) |> head()
#> [1]    0.000 1716.196 3290.379 4824.087 6745.759 7446.660

# 性能对比
bench::mark(
  original = cumulative_length(lines),
  optimized = cumulative_length_fast(lines)
)

2. terra 解决方案 #

利用terra的segments()拆分线段，结合distance()计算长度：

library(terra)
library(spData)

# 转换为SpatVector
lines_terra <- vect(seine[2, ])

cumulative_length_terra <- function(input) {
  # 拆分线段为单个片段
  segs <- segments(input)
  # 获取各片段长度
  seg_lengths <- distance(segs, along = TRUE)[, "length"]
  # 生成累积距离向量
  c(0, cumsum(seg_lengths))
}

# 验证结果
cumulative_length_terra(lines_terra) |> head()
#> [1]    0.000 1716.196 3290.379 4824.087 6745.759 7446.660

3. rsgeo 解决方案 #

rsgeo提供line_segment_lengths()直接获取线段各段长度，无需手动拆分：

library(rsgeo)
library(sf)
library(spData)

lines <- seine[2, ]

cumulative_length_rsgeo <- function(input) {
  # 转换为rsgeo Line对象
  line <- as_rsgeo(input)
  # 获取各线段长度
  seg_lengths <- line_segment_lengths(line)
  # 生成累积距离向量
  c(0, cumsum(seg_lengths))
}

# 验证结果
cumulative_length_rsgeo(lines) |> head()
#> [1]    0.000 1716.196 3290.379 4824.087 6745.759 7446.660

性能说明 #

上述三种方案均避免了原方法中循环创建大量几何对象的开销，其中sf+geos优化版和rsgeo的实现速度最快，terra方案也远优于原实现。

内容的提问来源于stack exchange，提问作者atsyplenkov