如何用R的sf包计算格网单元内道路网络的总长度？

核心思路 #

通过空间相交切割道路，再按格网单元分组统计切割后道路片段的总长度，具体步骤如下：

步骤1：加载依赖包并读取数据 #

先加载sf（空间处理）和dplyr（数据分组统计）：

library(sf)
library(dplyr)

读取道路（线要素）和格网（面要素）的shapefile：

# 读取道路数据（线）
roads <- st_read("path/to/your/roads.shp")
# 读取格网数据（面）
grid <- st_read("path/to/your/grid.shp")

步骤2：统一坐标系 #

确保两个空间图层的坐标系（CRS）完全一致，否则空间运算会出错：

# 检查当前CRS
st_crs(roads)
st_crs(grid)

# 如果不一致，将道路投影为格网的CRS
if (!st_crs(roads) == st_crs(grid)) {
  roads <- st_transform(roads, st_crs(grid))
}

步骤3：切割道路并计算长度 #

用st_intersection将道路切割为每个格网单元内的片段，然后计算每个片段的长度，再按格网的唯一标识分组求和：

# 若格网无唯一标识字段，先添加
grid <- grid %>% mutate(grid_id = row_number())

# 空间相交，得到每个格网内的道路片段
grid_road_intersect <- st_intersection(roads, grid)

# 计算每个道路片段的长度，并按格网ID分组求和
road_length_per_grid <- grid_road_intersect %>%
  mutate(segment_length = st_length(.)) %>%  # 计算单个片段长度
  as_tibble() %>%  # 转换为普通数据框简化统计
  group_by(grid_id) %>%
  summarise(total_road_length = sum(segment_length))

步骤4：合并结果到原始格网 #

将统计结果合并回原始格网，确保无道路的格网单元显示长度为0：

final_grid <- grid %>%
  left_join(road_length_per_grid, by = "grid_id") %>%
  mutate(total_road_length = replace_na(total_road_length, 0))

注意事项 #

• st_length()返回的长度单位由坐标系决定：投影坐标系（如UTM）返回米，地理坐标系（WGS84）返回度，可根据需求转换单位（比如用units::set_units()）
• 数据量较大时，建议先用st_make_valid()修复可能的无效几何：

roads <- st_make_valid(roads)
grid <- st_make_valid(grid)

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