求助：基于R的随机多边形模拟空间工作流实现方案

一、环境准备与数据加载 #

• 安装并加载核心依赖包：

install.packages(c("raster", "sf", "dplyr", "terra"))
library(raster)
library(sf)
library(dplyr)
library(terra)

• 加载目标数据并验证投影与面积：

# 加载Suit栅格（确保路径正确）
suit_rast <- raster("Suit.tif")
# 加载mpa多边形并统一投影
mpa_sf <- st_read("mpa.shp") %>% st_transform(crs = st_crs(suit_rast))

# 验证原mpa总面积是否为1465km²
total_mpa_area <- st_area(mpa_sf) %>% sum() %>% units::set_units(km²)
stopifnot(total_mpa_area == 1465)

二、随机多边形生成函数 #

生成1-3个总面积匹配、形状不规则且边缘锯齿状的多边形，采用随机点三角剖分实现自然不规则形态：

generate_random_polygons <- function(raster_extent, total_area_km2) {
  # 转换面积单位为平方米（适配32615投影的米单位）
  total_area_m2 <- total_area_km2 * 1e6
  # 随机确定本组多边形数量（1-3个）
  n_polys <- sample(1:3, 1)
  
  # 随机分配每个多边形的面积（确保总和匹配）
  area_ratios <- runif(n_polys)
  area_ratios <- area_ratios / sum(area_ratios)
  target_areas <- total_area_m2 * area_ratios
  
  # 逐个生成多边形
  poly_list <- lapply(target_areas, function(area) {
    # 生成随机点集（点数量控制锯齿程度）
    n_points <- sample(12:35, 1)
    points <- st_sample(raster_extent, size = n_points)
    # 构建Delaunay三角剖分并累加至目标面积
    tris <- st_triangulate(points)
    tri_areas <- st_area(tris)
    selected_tris <- c()
    current_area <- 0
    
    while(current_area < area * 0.95) {
      remaining <- setdiff(seq_along(tris), selected_tris)
      if(length(remaining) == 0) break
      pick <- sample(remaining, 1)
      selected_tris <- c(selected_tris, pick)
      current_area <- sum(st_area(tris[selected_tris]))
    }
    # 合并三角为单个多边形，微调面积偏差
    poly <- st_union(tris[selected_tris]) %>% st_cast("POLYGON")
    area_diff <- area - st_area(poly)
    if(abs(area_diff) > area * 0.05) {
      buffer_dist <- sqrt(abs(area_diff)/pi) * sign(area_diff)
      poly <- st_buffer(poly, dist = buffer_dist)
    }
    return(poly)
  })
  
  # 转换为sf对象返回
  do.call(rbind, lapply(poly_list, function(x) st_sf(geometry = x)))
}

三、批量模拟与指标计算 #

# 设置模拟次数
n_sim <- 100
sim_results <- list()

# 先计算原mpa内Suit>0.5的面积
mpa_suit_area <- suit_rast %>%
  mask(mpa_sf) %>%
  calc(function(x) ifelse(x > 0.5, 1, 0)) %>%
  cellStats(sum, na.rm = TRUE) %>%
  multiplyBy(res(suit_rast)[1] * res(suit_rast)[2]) %>%
  units::set_units(km²)

# 批量执行模拟
for(i in 1:n_sim) {
  # 生成随机多边形组
  rand_polys <- generate_random_polygons(st_as_sfc(st_bbox(suit_rast)), 1465)
  # 计算本组内每个多边形的Suit>0.5面积
  suit_areas <- lapply(1:nrow(rand_polys), function(j) {
    suit_rast %>%
      mask(rand_polys[j, ]) %>%
      calc(function(x) ifelse(x > 0.5, 1, 0)) %>%
      cellStats(sum, na.rm = TRUE) %>%
      multiplyBy(res(suit_rast)[1] * res(suit_rast)[2]) %>%
      units::set_units(km²)
  }) %>% unlist()
  
  # 存储结果
  sim_results[[i]] <- list(
    poly_count = nrow(rand_polys),
    total_suit_area = sum(suit_areas)
  )
  
  # 打印进度
  if(i %% 10 == 0) cat("完成模拟", i, "/", n_sim, "\n")
}

# 转换结果为数据框便于分析
sim_df <- do.call(rbind, lapply(sim_results, function(x) {
  data.frame(
    poly_count = x$poly_count,
    total_suit_area = as.numeric(x$total_suit_area)
  )
}))

四、选址合理性评估 #

通过对比原mpa与随机组的指标分布判断优劣：

# 计算原mpa在模拟分布中的分位数
mpa_quantile <- ecdf(sim_df$total_suit_area)(as.numeric(mpa_suit_area))

# 可视化对比
hist(sim_df$total_suit_area, main = "随机组Suit>0.5面积分布", 
     xlab = "面积(km²)", col = "lightblue")
abline(v = as.numeric(mpa_suit_area), col = "red", lwd = 2)
text(x = as.numeric(mpa_suit_area) + 40, y = max(hist(sim_df$total_suit_area, plot = FALSE)$counts), 
     labels = paste0("原mpa: ", round(as.numeric(mpa_suit_area), 1), "km²\n分位数: ", round(mpa_quantile, 3)), 
     col = "red")

• 若分位数>0.9：说明90%的随机组都达不到原mpa的目标面积，选址显著更优；
• 若分位数接近0.5：说明选址无明显优势，与随机选址效果相当。

注意事项 #

• 若计算速度慢，可替换raster为terra包，提升栅格处理效率；
• 可调整generate_random_polygons中的n_points范围，控制锯齿边缘的粗糙程度；
• 确保所有空间数据的CRS完全一致，避免投影偏差。

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