sf：生成具有最大距离条件的随机点答案

【问题标题】：sf: Generate random points with maximal distance conditionsf：生成具有最大距离条件的随机点
【发布时间】：2021-10-15 20:05:16
【问题描述】：

我想生成 100 个随机点，但在每个点周围使用大小为 1000 米的 st_buffer() 在点周围施加最大距离，并消除任何违规点。但是，在我的例子中：

library(sf)

# Data set creation
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")
st_bbox(df.laea)
# 

# Random simulation of 100 point inside df.laea extent
sim_study_area <- st_sample(st_as_sfc(st_bbox(df.laea)), 100) %>% # random points, as a list ...
  st_sf() 
border_area <- st_as_sfc(st_bbox(df.laea))%>% # random points, as a list ...
  st_sf() 

# I'd like to imposed a maximal distance of 1000 meters around points and for this:

i <- 1 # iterator start

buffer_size <- 1000 # minimal distance to be enforced (in meters)

repeat( {
    
  #  create buffer around i-th point
  buffer <- st_buffer(sim_study_area[i,], buffer_size) 
  
  offending <- sim_study_area %>%  # start with the intersection of master points... 
    st_intersects(buffer, sparse = F) # ... and the buffer, as a vector
  
  # i-th point is not really offending 
  offending[i] <- TRUE
  
  # if there are any offending points left - re-assign the master points
  sim_study_area <- sim_study_area[offending,] 
  
  if ( i >= nrow(sim_study_area)) {
      # the end was reached; no more points to process
      break 
    } else {
      # rinse & repeat
      i <- i + 1 
    }
  
} )

# Visualizantion of points create with the offending condition:
simulation_area <- ggplot() +
  geom_sf(data = border_area, col = 'gray40', fill = NA, lwd = 1) +
  geom_sf(data = sim_study_area, pch = 3, col = 'red', alpha = 0.67) +
  theme_bw() 
plot(simulation_area)

结果不好，因为没有 100 分，我不知道如何解决。

请问有什么想法吗？

提前致谢，

亚历山大

【问题讨论】：

很难确定满足特定标准的随机点的确切数量 - 点的随机性使它们有点不可预测。我怀疑这必须通过反复试验来完成。

标签： r sf sp

【解决方案1】：

我认为最简单的解决方案是采用R包spatstat中定义的采样函数之一。例如：

# packages
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

# create data
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"
)

现在我们使用简单的顺序抑制过程进行采样。更多详情请查看?spatstat.core::rSSI。

sampled_points <- st_sample(
  x = st_as_sfc(st_bbox(df.laea)),
  type = "SSI",
  r = 1000, # threshold distance (in metres)
  n = 100 # number of points
)

# Check result
par(mar = rep(0, 4))
plot(st_as_sfc(st_bbox(df.laea)), reset = FALSE)
plot(sampled_points, add = TRUE, pch = 16)

# Estimate all distances
all_distances <- st_distance(sampled_points)
all_distances[1:5, 1:5]
#>           [,1]      [,2]      [,3]      [,4]      [,5]
#> [1,]      0.00  57735.67 183205.74 189381.50  81079.79
#> [2,]  57735.67      0.00 153892.93 143755.73  61475.85
#> [3,] 183205.74 153892.93      0.00  62696.68 213379.39
#> [4,] 189381.50 143755.73  62696.68      0.00 194237.12
#> [5,]  81079.79  61475.85 213379.39 194237.12      0.00

# Check they are all greater than 1000
sum(all_distances < 1000)
#> [1] 100 # since the diagonal is full of 100 zeros

^{由reprex package (v2.0.0) 于 2021-08-12 创建}

查看here（特别是Baddeley 教授的回答）、其中的参考资料以及st_sample 的帮助页面了解更多详情。

【讨论】：