如何通过与栅格(raster)数据相交生成sf格式空间矢量数据

针对矢量街道数据与栅格坡度数据的交叉分析需求（生成带坡度列的矢量街道数据，找出最陡路段），无需使用st_intersection（该函数适用于矢量间交互），正确的做法是通过栅格值提取将坡度信息关联到街道要素上。以下是完整可复现的实现方案：

完整代码 #

library(elevatr) 
library(terra)
library(geobr)
library(osmdata)
library(sf)
library(dplyr)

# 获取市政边界矢量
municipality_shape <- read_municipality(code_muni = 3305802)

# 获取地形栅格并转换为terra格式
elev_raster <- get_elev_raster(locations = municipality_shape, 
                               z = 10, prj = "EPSG:4674") 
obj_raster <- rast(elev_raster) 

# 计算坡度栅格（百分比单位更直观，也可保留弧度）
slope_raster <- terrain(obj_raster, "slope", unit = "percent")

# 获取主干道矢量数据
big_streets <- getbb("Teresópolis") |>
  opq() |>
  add_osm_feature(key = "highway", value = "primary") |>
  osmdata_sf()

# 整理街道为sf对象，统一坐标系
streets_sf <- big_streets$osm_lines
streets_sf <- st_transform(streets_sf, crs = crs(slope_raster))

# 提取每条街道覆盖范围内的最大坡度（匹配"最陡路段"需求）
slope_results <- terra::extract(slope_raster, vect(streets_sf), fun = max, na.rm = TRUE)

# 合并坡度列到街道数据
streets_with_slope <- cbind(streets_sf, max_slope = slope_results[,2])

# 筛选坡度最陡的路段（示例：取坡度前10%的路段）
steepest_streets <- streets_with_slope |>
  arrange(desc(max_slope)) |>
  slice_head(prop = 0.1)

# 查看结果
head(streets_with_slope[,c("name", "max_slope", "geometry")])

关键细节说明 #

• 坐标系统一：矢量与栅格数据的坐标系必须完全一致，否则提取会出现错误。使用st_transform()将街道数据转换为与坡度栅格相同的坐标系。
• 栅格提取逻辑：terra::extract()支持对线要素提取栅格统计值：
  • 指定fun = max可获取每条街道路段覆盖范围内的最大坡度值，精准匹配"找出最陡路段"的需求；
  • 若需要路段的平均坡度，可将fun改为mean。
• 结果输出：合并后的streets_with_slope是带max_slope列的sf矢量数据，可直接用于后续的筛选、可视化或进一步分析。

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