设置锚点并从Shapefile创建多边形：适配代码遇几何有效性报错

错误原因分析 #

1. CRS不匹配警告：你使用的锚点坐标是投影坐标系数值，但the_lake可能是经纬度坐标系，导致系统误判坐标范围，触发警告。
2. 几何无效错误：旋转后的条带多边形或原湖泊几何存在自相交、重复顶点问题，导致st_intersection无法正常执行。

修正方案与代码 #

1. 统一坐标系并修复几何有效性 #

先确保所有空间数据使用同一投影坐标系，同时修复几何图形的有效性：

# 加载库
library(sf)
library(dplyr)

# 读取数据并统一CRS
the_lake <- readRDS('vallejo.RDS')
# 转换为适配Vallejo的UTM投影（EPSG:32610），修复湖泊几何
target_crs <- 32610
the_lake <- the_lake %>% 
  st_transform(target_crs) %>% 
  st_make_valid()

# 将原经纬度锚点转换为目标投影坐标
anchor_latlon <- st_sfc(st_point(c(-122.244700, 38.067690)), crs = 4326)
the_anchor <- st_transform(anchor_latlon, target_crs) %>% st_coordinates() %>% as.vector()

# 条带生成函数：调整默认CRS为目标投影
get_strips <- function(n = 10, 
                       strip_width = 1e3, 
                       strip_length = 1e5,
                       crs = target_crs
                       ){
  baseline <- matrix(c(-strip_length/2, strip_length/2, 0, 0), ncol = 2) %>%
    st_linestring()
  basestrip  <- baseline %>% 
    st_buffer(strip_width/2, endCapStyle = 'FLAT')
  st_sf(strip_id = sprintf('strip_%02.f', 1:n),
        geometry = Map(1:n,
                       f = \(index) basestrip - c(0, index * strip_width - strip_width/2)
                       ) %>%
          st_sfc() %>% st_cast('MULTIPOLYGON') %>% st_set_crs(crs)
        )
}

# 旋转函数：新增旋转后几何修复
rotate_feature <- function(feature,
                           angle = 0,
                           anchor = c(0, 0)
                           ){
  crs <- st_crs(feature)
  a <- -pi * angle/180
  rotmat <- matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
  st_geometry(feature) <- (st_geometry(feature) - st_point(c(0,0))) * rotmat + st_point(anchor)
  feature %>% st_set_crs(crs) %>% st_make_valid()
}

# 生成条带并裁剪
the_strips <- 
  get_strips(crs = target_crs) %>%
  rotate_feature(angle = 10, anchor = the_anchor)

the_cropped_strips <- the_strips %>% 
  st_make_valid() %>%
  st_intersection(the_lake)

# 提取相交角点：明确转换为点类型
the_corners <- 
  st_intersection(
    st_cast(the_strips %>% st_make_valid(), 'MULTILINESTRING'),
    st_cast(the_lake, 'MULTILINESTRING')
  ) %>%
  st_cast('POINT')

关键修正点 #

• 统一使用UTM投影（EPSG:32610），避免经纬度与投影坐标混用的问题。
• 新增st_make_valid()调用，修复所有空间几何的自相交、重复顶点问题。
• 将锚点从经纬度转换为目标投影坐标，确保旋转逻辑正确。
• 提取角点时明确转换为POINT类型，避免混合几何类型问题。

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