如何基于bay.RDS制作指定位置的50米宽跨岸多边形并提取坐标？

1. 加载依赖包与数据 #

首先加载空间分析所需的包，再加载本地的bay.RDS文件：

library(tmap)
library(leaflet)
library(mapview)
library(sf) # 补充空间操作核心包

# 加载本地的bay.RDS文件（需确保文件在当前工作目录）
bay <- readRDS('bay.RDS')
mapview(bay) # 预览底图

2. 定义多个相邻目标点位 #

指定需要生成多边形的相邻点位（可根据需求调整坐标数量与位置）：

# 定义3个相邻点位示例
pts <- st_sfc(
  st_point(c(-122.91, 38.17)),
  st_point(c(-122.908, 38.17)),
  st_point(c(-122.906, 38.17)),
  crs = 4326
)

# 预览点位与底图叠加效果
mapview(bay) + mapview(pts, color = 'red')

3. 批量生成岸到岸的50米宽多边形 #

核心逻辑：先转平面坐标系（米为单位），基于点位生成垂直于岸线的宽缓冲区，再裁剪到湾区范围得到岸到岸的独立多边形：

# 转换为UTM平面坐标系（适配旧金山湾区的EPSG:32610，单位为米）
bay_utm <- st_transform(bay, crs = 32610)
pts_utm <- st_transform(pts, crs = 32610)

# 提取湾区的边界线
bay_edge <- st_cast(bay_utm, "LINESTRING")

# 定义单点位生成多边形的函数
create_cross_poly <- function(pt, edge, width = 50) {
  # 计算点到岸线的垂直方向
  nearest_pt <- st_nearest_points(pt, edge)[[1]][2]
  direction <- st_coordinates(pt - nearest_pt)[1,]
  
  # 生成横跨两岸的线段（延伸足够长度确保覆盖湾区）
  cross_line <- st_sfc(st_linestring(rbind(
    st_coordinates(pt) + direction * 10000,
    st_coordinates(pt) - direction * 10000
  )), crs = st_crs(edge))
  
  # 生成50米宽缓冲区并裁剪到湾区范围
  buffer <- st_buffer(cross_line, width / 2)
  poly <- st_intersection(buffer, bay_utm)
  return(st_cast(poly, "POLYGON"))
}

# 批量生成所有点位对应的多边形
cross_polys <- lapply(pts_utm, function(x) create_cross_poly(x, bay_edge))
cross_polys <- st_sf(geom = do.call(c, cross_polys))

# 添加用于颜色编码的属性列（替换为你的实际业务数据）
cross_polys$score <- c(15, 25, 35)

# 转换回WGS84经纬度坐标系（EPSG:4326）
cross_polys_wgs84 <- st_transform(cross_polys, crs = 4326)

# 预览效果：按score列颜色编码多边形
mapview(bay) + mapview(cross_polys_wgs84, zcol = "score")

4. 提取多边形坐标并合并数据集 #

将多边形的经纬度坐标提取为表格格式，可直接合并到你的业务数据中：

# 定义提取单个多边形坐标的函数
get_poly_coords <- function(poly_row) {
  coords <- st_coordinates(poly_row$geom)
  return(data.frame(
    poly_id = rep(poly_row$row.names, nrow(coords)),
    lon = coords[,1],
    lat = coords[,2],
    score = rep(poly_row$score, nrow(coords))
  ))
}

# 批量提取所有多边形的坐标
coords_df <- do.call(rbind, lapply(1:nrow(cross_polys_wgs84), function(i) {
  get_poly_coords(cross_polys_wgs84[i,])
}))

# 查看提取结果
head(coords_df)

# 合并到你的原始数据集（假设原始数据集为raw_df，包含对应poly_id）
# merged_data <- merge(raw_df, coords_df, by = "poly_id")

效果说明 #

生成的多边形为相邻但独立的岸到岸区域，可通过自定义属性列（如示例中的score）实现颜色编码，提取的坐标为标准经纬度格式，可直接用于后续分析或可视化。

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