R sf包：学校点450m缓冲区重叠时对应250m缓冲区合并方法

现有地图上的点要素为学校点位，每个学校点位对应生成两个缓冲区：250m缓冲区和450m缓冲区。需求为当点位的450m缓冲区存在重叠时，将这些点位归为同一个统计单元，同时保留各单元覆盖的几何范围。如果直接对250m缓冲区和450m缓冲区分别调用st_union函数合并，会因450m缓冲区的重叠范围更广，导致两类缓冲区合并后得到的单元数量不一致，无法匹配。

核心需求为：只要学校的450m缓冲区重叠合并，对应的250m缓冲区也同步合并为同一个单元，最终统计商店分别落在学校250m缓冲区、450m缓冲区内的次数。

核心逻辑是先以450m缓冲区的重叠关系为基准给学校分组，再按组合并两类缓冲区，从根源上保证两类缓冲区合并后单元数量完全匹配。

实现步骤 #

• 给每个学校点位生成唯一ID，用于后续分组关联
• 基于450m缓冲区的邻接关系计算连通分量，每个连通分量对应一个统计单元，所有属于同一个连通分量的学校归为一组
• 按分组ID分别合并250m和450m缓冲区，得到的两个结果集单元数量完全一致，一一对应

完整代码 #

library(tidycensus)
library(sf)
library(tmap)
library(rio)
library(dplyr)
library(igraph)

# 基础底图数据读取
county.sf <- get_acs(state = "MO",
                     county = c( "St. Louis City"),
                     geography = "tract",
                     variables = "B03002_001", 
                     output="wide", 
                     geometry = TRUE) %>%
  sf::st_transform(crs = "ESRI:102003")

# 学校数据读取与坐标系转换
school <- read.csv("C:\\myfile1.csv")
school.sf <- st_as_sf(school, coords = c("long", "lat"), crs = "epsg:4326") 
school.sf.utm <- st_transform(school.sf, crs = "ESRI:102003")

# 商店数据读取与坐标系转换
store <- import("C:\\myfile2.csv")
store.sf <- st_as_sf(store, coords = c("XCoord", "YCoord"), crs = "ESRI:102696") 
store.sf.utm <- st_transform(store.sf, crs = "ESRI:102003")

# 生成两类缓冲区
elem.buff <- st_buffer(school.sf.utm, 250)     
elem.buff2 <- st_buffer(school.sf.utm, 450) 

# 按450m缓冲区重叠关系生成分组ID
elem.buff2$school_id <- 1:nrow(elem.buff2)
buff2_intersect <- st_intersects(elem.buff2, elem.buff2)
group_ids <- components(graph_from_adj_list(buff2_intersect))$membership
elem.buff$group_id <- group_ids
elem.buff2$group_id <- group_ids

# 按分组合并两类缓冲区，保证单元数完全匹配
merged_buff250 <- elem.buff %>% 
  group_by(group_id) %>% 
  summarise(geometry = st_union(geometry), .groups = "drop")
merged_buff450 <- elem.buff2 %>% 
  group_by(group_id) %>% 
  summarise(geometry = st_union(geometry), .groups = "drop")

# 统计每个单元内的商店数量
merged_buff250$store_count_250m <- lengths(st_intersects(merged_buff250, store.sf.utm))
merged_buff450$store_count_450m <- lengths(st_intersects(merged_buff450, store.sf.utm))

# 合并结果可视化
ex.map <- tm_shape(county.sf) +
  tm_polygons() + 
  tm_shape(merged_buff250) +
  tm_borders(col="red") +  
  tm_shape(merged_buff450) +
  tm_borders(col="blue") + 
  tm_shape(school.sf.utm) +
  tm_dots(col = "red") +
  tm_shape(store.sf.utm) +
  tm_dots() 
print(ex.map)

结果说明 #

合并后的merged_buff250和merged_buff450行数完全相同，同一个group_id对应同一个统计单元，只要450m缓冲区重叠的学校就会被分到同一组，对应的250m缓冲区也会合并到同一单元，完全匹配需求。最终统计得到的store_count_250m和store_count_450m分别对应每个单元内两类缓冲区覆盖的商店数量。

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