如何在R中实现线串(linestring)到多边形(polygon)的完全贴合？

我有两个SHP文件：一个包含多条线串（铁路数据），另一个包含多个多边形（美国县界数据），两者完全没有交集。我的目标是将线串贴合到最近的多边形边界，使其产生交集，进而计算每条线串在各多边形内的长度。

我尝试用sf包的st_snap函数实现，但未得到预期结果，于是用简单数据测试该函数：

测试代码与初始状态 #

library(dplyr)
library(tidyverse)
library(sf)
library(units)
library(ggplot2)

# 创建测试多边形与线串
poly = st_polygon(list(cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))))
lines = st_multilinestring(list(
 cbind(c(0, 1), c(1.01, 1.05)),
 cbind(c(0, 1), c(-0.01, -.05)),
 cbind(c(1.025, 1.05, 1.025), c(1.05, .5, -.05))
))

# 验证初始无交集
t <- st_intersection(poly, lines)
st_is_empty(t)
# [1] TRUE

# 可视化初始状态
ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)

![线（绿色）&多边形（红色）]

我的预期是线串贴合到多边形最近的边界上，与多边形产生交集。

st_snap的实际表现 #

当tolerance=0.5时 #

snapped1 = st_snap(lines, poly, tolerance=0.5)
print(snapped1)
# MULTILINESTRING ((0 1, 1 1), (0 0, 1 0), (1 1, 1.05 0.5, 1 0))

ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)+
    geom_sf(data = snapped1, col = "blue", alpha = 0.5)

![蓝色贴合线效果]

当tolerance=1.005时 #

snapped2 = st_snap(lines, poly, tolerance=1.005)
print(snapped2)
# MULTILINESTRING ((1 1, 0 1, 0 0, 1 0), (0 1, 1 1, 1 0, 0 0), (0 1, 1 1, 1.05 0.5, 1 0, 0 0))

ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)+
    geom_sf(data = snapped2, col = "blue", alpha = 0.5)

结果包含预期的(0 1, 1 1, 1 0, 0 0)，但同时存在多余线段。
![蓝色贴合线效果]

当tolerance=1.1时 #

snapped3 = st_snap(lines, poly, tolerance=1.1)
print(snapped3)
# MULTILINESTRING ((1 1, 0 1, 0 0, 1 0), (0 1, 1 1, 1 0, 0 0), (1 1, 1 0, 0 0, 0 1))

ggplot()+
    geom_sf(data = lines, col = "green")+
    geom_sf(data = poly, col = "red", fill = NA)+
    geom_sf(data = snapped3, col = "blue", alpha = 0.5)

结果同样包含预期部分，但存在多余线段，更大容差会得到相同结果。

st_snap工作原理解析 #

st_snap的核心逻辑是将目标几何的顶点捕捉到参考几何的顶点或边上，但它会保留原几何的拓扑结构：

• 当容差小于线到多边形的距离时，顶点不会被捕捉，结果无变化；
• 当容差刚好覆盖线到多边形的距离时，顶点会被捕捉到最近的边界点；
• 当容差过大时，线的顶点会被捕捉到多边形的多个顶点，甚至会连接不同的边界点，生成绕多边形的多余线段——这就是你遇到的问题。

实现预期结果的方法 #

方法一：st_snap结合裁剪（批量处理友好） #

先计算线到多边形的最小距离，设置刚好覆盖该距离的容差，避免过度捕捉，再用多边形裁剪掉多余部分：

# 计算每条线到多边形的最小距离
line_distances <- st_distance(lines, poly)
# 设置容差为最大距离+微小值，确保刚好能捕捉
tolerance_val <- max(line_distances) + 0.001

# 执行捕捉
snapped_lines <- st_snap(lines, poly, tolerance = tolerance_val)
# 裁剪到多边形内部
clipped_lines <- st_intersection(snapped_lines, poly)

# 查看结果
print(clipped_lines)
ggplot() +
  geom_sf(data = lines, col = "green") +
  geom_sf(data = poly, col = "red", fill = NA) +
  geom_sf(data = clipped_lines, col = "blue", alpha = 0.5)

方法二：手动投影线串到最近边界（精准控制） #

将多边形转为线，找到每条线到多边形的最近边，手动投影线的端点到该边上，生成新线串后再裁剪：

# 将多边形转为线
poly_line <- st_cast(poly, "LINESTRING")

# 定义投影函数
project_line_to_poly <- function(line, poly_line) {
  # 找到线与多边形边的最近点对
  nearest_points <- st_nearest_points(line, poly_line)
  # 提取投影后的端点（取每个点对的第二个点，即多边形边上的点）
  proj_start <- st_cast(nearest_points[1,], "POINT")[2]
  proj_end <- st_cast(nearest_points[2,], "POINT")[2]
  # 生成新的线串
  st_linestring(rbind(st_coordinates(proj_start), st_coordinates(proj_end)))
}

# 拆分多线串为单个线串，逐个投影
lines_single <- st_cast(lines, "LINESTRING")
projected_lines <- lapply(lines_single, project_line_to_poly, poly_line = poly_line)
projected_lines <- st_sfc(projected_lines, crs = st_crs(lines))

# 裁剪到多边形内部
clipped_projected <- st_intersection(projected_lines, poly)

# 查看结果
print(clipped_projected)
ggplot() +
  geom_sf(data = lines, col = "green") +
  geom_sf(data = poly, col = "red", fill = NA) +
  geom_sf(data = clipped_projected, col = "blue", alpha = 0.5)

后续计算长度 #

得到贴合后的线串后，可以用st_length计算每条线在多边形内的长度，结合分组统计：

# 假设已将线串与多边形关联（例如通过st_join）
length_stats <- clipped_projected %>%
  st_join(st_as_sf(poly)) %>%
  mutate(length = st_length(.)) %>%
  group_by(your_poly_id_column) %>%
  summarise(total_length = sum(length))

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