如何使用R从道路shapefile生成1×1km道路密度栅格

前置依赖 #

需要用到三个R空间处理包：sf（矢量数据处理）、terra（栅格数据处理）、exactextractr（栅格与矢量精确交叠统计），相比传统矢量叠加的方案，该组合处理大区域数据效率更高。

步骤1：数据读取与投影校正 #

# 首次使用先安装依赖包
install.packages(c("sf", "terra", "exactextractr"))
# 加载包
library(sf)
library(terra)
library(exactextractr)

# 读取道路shapefile，替换为你本地的文件路径
road_sf <- st_read("your_road_file.shp")
# 转换为投影单位为米的等面积投影，此处用适配美国全境的US National Atlas Equal Area投影
# 也可替换为你常用的美国东部区域投影，只要单位为米即可
road_sf_proj <- st_transform(road_sf, crs = "EPSG:2163")

注意：必须使用单位为米的平面投影，否则1km格网大小、长度统计都会出现偏差，不能直接用WGS84等经纬度投影计算。

步骤2：生成1km分辨率栅格模板 #

# 获取道路数据的边界范围
road_bbox <- st_bbox(road_sf_proj)
# 生成空栅格，res参数设置为1000即对应1×1km分辨率
raster_temp <- rast(road_bbox, res = c(1000, 1000), crs = st_crs(road_sf_proj)$wkt)
# 给每个栅格赋唯一ID用于后续统计
values(raster_temp) <- 1:ncell(raster_temp)

步骤3：单像元道路长度统计 #

因为你的shapefile已经自带道路长度属性列，可直接调用该字段统计，无需重复计算长度：

# 假设道路长度属性列名为length，替换为你实际的列名
# 统计每个像元内所有道路的长度总和
density_stat <- exact_extract(
  x = raster_temp,
  y = road_sf_proj,
  fun = "sum",
  weights = "length",
  summarize_df = TRUE
)
# 把统计结果赋值到栅格模板
values(raster_temp) <- density_stat$sum
# 无道路的像元值替换为0
raster_temp[is.na(raster_temp)] <- 0
# 导出结果栅格，像元值单位为米/平方千米，除以1000即可转换为千米/平方千米
writeRaster(raster_temp, "us_east_road_density_1km.tif", overwrite = TRUE)

若数据量较小，也可以用矢量叠加方案：先生成1km矢量格网，用st_intersection裁剪道路与格网，再按格网ID分组求和长度即可，该方案逻辑更直观但大区域数据下运行速度较慢。

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