如何在R中将指定Shapefile随机划分为n个不重叠子区域？

以下是三种在R中实现随机分区并导出Shapefile的方法，基于你已有的代码基础：

方法一：基于Voronoi图（泰森多边形）的随机分区 #

通过在目标区域内生成随机点，构建Voronoi多边形后与原区域裁剪，得到不重叠的子区域。

代码实现 #

library(sf)
library(dplyr)
library(tmap) # 可选，用于可视化

# 基础代码：加载并处理目标区域
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>% 
  st_transform(st_crs(4326)) %>% 
  st_cast('POLYGON') %>%
  filter(NAME %in% c("Cleveland", "Gaston", "Rutherford", "Burke")) %>%
  st_union() %>% # 合并为单个多边形
  st_sf() # 转换为sf对象

n <- 8 # 子区域数量

# 1. 在目标区域内生成n个随机点
set.seed(123) # 设置随机种子保证可复现
random_points <- st_sample(nc, size = n, type = "random") %>%
  st_sf() %>%
  mutate(id = 1:n)

# 2. 生成Voronoi图并裁剪到原区域边界
voronoi <- st_voronoi(st_union(random_points)) %>%
  st_collection_extract("POLYGON") %>%
  st_sf() %>%
  mutate(id = 1:n) %>%
  st_intersection(nc) # 裁剪到原区域

# 3. 可视化结果（可选）
tm_shape(voronoi) + tm_polygons(col = "id", alpha = 0.7) + tm_borders(col = "black")

# 4. 导出为Shapefile
st_write(voronoi, "random_voronoi_subregions.shp", delete_layer = TRUE)

说明 #

• 优点：子区域边界自然，适合需要随机分散的分区场景
• 注意：如果随机点分布不均，子区域面积差异可能较大，可通过分层随机采样优化

方法二：使用spatialEco包生成随机多边形分区 #

spatialEco包的randomPolygons函数可直接在目标区域内生成指定数量的不重叠随机多边形。

代码实现 #

library(sf)
library(dplyr)
library(spatialEco)

# 基础代码同上
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>% 
  st_transform(st_crs(4326)) %>% 
  st_cast('POLYGON') %>%
  filter(NAME %in% c("Cleveland", "Gaston", "Rutherford", "Burke")) %>%
  st_union() %>%
  st_sf()

n <- 8

# 转换为sp对象（spatialEco依赖sp格式）
nc_sp <- as(nc, "Spatial")

# 生成随机不重叠多边形
set.seed(123)
random_polys <- randomPolygons(nc_sp, n = n, min = 0.1, max = 0.5, 
                               sp = TRUE, prob = NULL) %>%
  st_as_sf() %>% # 转换回sf格式
  mutate(subregion_id = 1:n)

# 可视化（可选）
plot(random_polys["subregion_id"], border = "black")

# 导出Shapefile
st_write(random_polys, "random_polygons_subregions.shp", delete_layer = TRUE)

说明 #

• min和max参数控制子区域面积的最小/最大比例（相对于原区域），可按需调整
• 优点：直接生成符合面积要求的随机多边形，无需额外裁剪

方法三：基于规则网格的随机合并分区 #

先将区域划分为细粒度规则网格，再随机合并网格单元得到指定数量的子区域，适合需要面积相对均匀的场景。

代码实现 #

library(sf)
library(dplyr)
library(tmap)

# 基础代码同上
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>% 
  st_transform(st_crs(4326)) %>% 
  st_cast('POLYGON') %>%
  filter(NAME %in% c("Cleveland", "Gaston", "Rutherford", "Burke")) %>%
  st_union() %>%
  st_sf()

n <- 8
grid_size <- 0.02 # 网格单元大小，根据区域尺度调整

# 1. 创建规则网格
grid <- st_make_grid(nc, cellsize = grid_size, what = "polygons") %>%
  st_sf() %>%
  st_intersection(nc) %>%
  mutate(grid_id = row_number())

# 2. 随机分配网格单元到n个子区域
set.seed(123)
grid <- grid %>%
  mutate(subregion_id = sample(1:n, nrow(grid), replace = TRUE))

# 3. 合并同一子区域的网格单元
merged_subregions <- grid %>%
  group_by(subregion_id) %>%
  summarise(geometry = st_union(geometry)) %>%
  st_cast("POLYGON") # 确保是单一多边形

# 4. 可视化
tm_shape(merged_subregions) + tm_polygons(col = "subregion_id", alpha = 0.7) + tm_borders(col = "black")

# 5. 导出Shapefile
st_write(merged_subregions, "merged_grid_subregions.shp", delete_layer = TRUE)

说明 #

• 调整grid_size可以控制网格精细度，网格越细，合并后的子区域边界越自然
• 优点：子区域面积相对均匀，适合需要大致均衡分区的场景

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