使用R的sf包生成满足1000米间距要求的100个随机点问题求助

前置检查 #

首先确认你的研究区面积是否足以放下100个满足最小间距要求的点：按六边形最密铺规则计算，100个间隔1000米的点至少需要约8.66e7平方米（即86.6平方公里）的区域。可运行st_area(st_as_sfc(st_bbox(df.laea)))查看研究区面积，如果小于阈值，不可能生成足够的点，需要先扩大研究区范围。

方案1：修改原有逻辑，补足采样量 #

原有代码仅采样100个初始点，剔除不符合要求的点后自然达不到目标数量，调整为初始采样多倍目标数量的点，不足时自动补充采样：

library(sf)
library(dplyr)

# 原有数据处理逻辑保留
set.seed(1)
df <- data.frame(
  gr = c(rep("a",5),rep("b",5)),
  x  = rnorm(10),
  y = rnorm(10)
)
df <- st_as_sf(df,coords = c("x","y"),remove = F, crs = 4326)
df.laea = st_transform(df, 
                       crs = "+proj=laea +x_0=4600000 +y_0=4600000 +lon_0=0.13 +lat_0=0.24 +datum=WGS84 +units=m")
study_bbox <- st_as_sfc(st_bbox(df.laea))
buffer_size <- 1000 
target_n <- 100

# 初始采样5倍目标数量的点，提高初始合格概率
init_points <- st_sample(study_bbox, target_n * 5) %>% st_sf()
valid_points <- st_sf(geometry = st_sfc(crs = st_crs(df.laea)))

i <- 1
while(nrow(valid_points) < target_n && i <= nrow(init_points)){
  current_pt <- init_points[i,]
  # 检查当前点和已有合格点的距离是否都大于等于阈值
  if(nrow(valid_points) == 0 || min(st_distance(current_pt, valid_points)) >= buffer_size){
    valid_points <- rbind(valid_points, current_pt)
  }
  i <- i + 1
  # 如果初始点用完还没凑够，补充新的采样点
  if(i > nrow(init_points) && nrow(valid_points) < target_n){
    add_points <- st_sample(study_bbox, target_n * 2) %>% st_sf()
    init_points <- rbind(init_points, add_points)
  }
}

# 可视化验证
library(ggplot2)
ggplot() +
  geom_sf(data = study_bbox %>% st_sf(), col = 'gray40', fill = NA, lwd = 1) +
  geom_sf(data = valid_points, pch = 3, col = 'red', alpha = 0.67) +
  theme_bw()

方案2：使用专用采样函数（效率更高） #

直接用spatstat包的rSSI函数，原生支持生成满足最小距离要求的随机点，无需手动写迭代逻辑：

library(sf)
library(spatstat)

# 原有数据处理逻辑
set.seed(1)
df <- data.frame(
  gr = c(rep("a",5),rep("b",5)),
  x  = rnorm(10),
  y = rnorm(10)
)
df <- st_as_sf(df,coords = c("x","y"),remove = F, crs = 4326)
df.laea = st_transform(df, 
                       crs = "+proj=laea +x_0=4600000 +y_0=4600000 +lon_0=0.13 +lat_0=0.24 +datum=WGS84 +units=m")
study_bbox <- st_bbox(df.laea)
buffer_size <- 1000 
target_n <- 100

# 转成spatstat支持的窗口对象
study_owin <- owin(xrange = c(study_bbox$xmin, study_bbox$xmax),
                   yrange = c(study_bbox$ymin, study_bbox$ymax))

# 生成满足最小距离的随机点，giveup为最大尝试次数，调大可提高凑够点数的概率
ss_points <- rSSI(r = buffer_size, n = target_n, win = study_owin, giveup = 10000)

# 转成sf对象
valid_points <- st_as_sf(as(ss_points, "SpatialPoints")) %>% 
  st_set_crs(st_crs(df.laea))

# 可视化
library(ggplot2)
ggplot() +
  geom_sf(data = st_as_sfc(study_bbox) %>% st_sf(), col = 'gray40', fill = NA, lwd = 1) +
  geom_sf(data = valid_points, pch = 3, col = 'red', alpha = 0.67) +
  theme_bw()

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