R中沿县界生成等距垂线技术求助：排除州界与共享边界去重

以下是针对你需求的分步实现方案，核心是先精准筛选内部共享边界，再进行等距采样和垂线计算：

1. 筛选内部共享边界（排除州外边界） #

首先区分仅属于单个县的州外边界和被两个县共享的内部边界：

library(sf)
library(dplyr)

# 假设你的县数据框名为counties_sf，包含唯一标识列county_id
# 提取每个县的边界并拆分为单个LINESTRING
counties_lines <- st_cast(st_boundary(counties_sf), "LINESTRING") %>%
  mutate(county_id = counties_sf$county_id)

# 统一线段方向（消除共享边界的方向差异），生成唯一标识字符串
counties_lines <- counties_lines %>%
  mutate(norm_geom = st_normalize(geometry),
         geom_wkt = st_as_text(norm_geom))

# 筛选出被两个县共享的边界（出现次数为2），并去重保留单条
shared_internal_lines <- counties_lines %>%
  group_by(geom_wkt) %>%
  filter(n() == 2) %>%
  ungroup() %>%
  distinct(geom_wkt, .keep_all = TRUE)

2. 等距采样共享边界 #

用st_line_sample可以精准控制采样间隔（需确保数据为米制投影坐标系，如UTM；若为地理坐标系先转换）：

# 转换为米制投影（示例用UTM 10N，EPSG:32610，根据你的州位置调整）
if (st_is_longlat(shared_internal_lines)) {
  shared_internal_lines <- st_transform(shared_internal_lines, crs = 32610)
}

# 按1000米间隔采样（density为每米采样数，1/1000即每1000米1个点）
sample_points <- st_line_sample(shared_internal_lines$norm_geom, density = 1/1000) %>%
  st_cast("POINT") %>%
  st_sf() %>%
  # 关联对应的原始边界线段
  mutate(line_idx = rep(1:nrow(shared_internal_lines), each = 
                          st_length(shared_internal_lines$norm_geom) %>%
                          as.numeric() %>%
                          floor() %>%
                          divide_by(1000) %>%
                          add(1)))

3. 计算采样点的垂线方向 #

通过计算边界线段的切线方向，旋转90度得到垂线方向（以下示例取左侧垂线方向）：

# 定义函数：输入线段和点，返回垂线方向的单位向量
get_perpendicular <- function(line, point) {
  coords <- st_coordinates(line)
  dir_vec <- coords[nrow(coords), 1:2] - coords[1, 1:2]
  dir_vec <- dir_vec / sqrt(sum(dir_vec^2)) # 单位化方向向量
  # 旋转90度得到左侧垂线，取反则为右侧
  c(-dir_vec[2], dir_vec[1])
}

# 批量计算每个采样点的垂线方向
sample_points <- sample_points %>%
  rowwise() %>%
  mutate(perp_x = get_perpendicular(shared_internal_lines$norm_geom[line_idx], geometry)[1],
         perp_y = get_perpendicular(shared_internal_lines$norm_geom[line_idx], geometry)[2]) %>%
  ungroup()

解决st_segmentize报错的小技巧 #

如果必须用st_segmentize，先处理无效几何并拆分线段：

# 修复无效几何并确保是单个LINESTRING
valid_lines <- shared_internal_lines %>%
  st_make_valid() %>%
  st_cast("LINESTRING")

# 指定最大段长（单位与投影一致，此处为米）
segmentized_lines <- st_segmentize(valid_lines$norm_geom, max_segment_length = 1000)

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