sf: Generate random points with maximal distance condition

Question

I'd like to generate 100 random points but imposed a maximal distance around points using st_buffer() of size 1000 meters around each point, and eliminating any offending points. But, in my example:

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)

It's not OK result because a don't have 100 points and I don't know how I can fix it.

Please any ideas?

Thanks in advance,

Alexandre

Answer 1

I think that the easiest solution is to adopt one of the sampling functions defined in the R package spatstat. For example:

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

Now we sample with a Simple Sequential Inhibition Process. Check ?spatstat.core::rSSI for more details.

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

^{Created on 2021-08-12 by the reprex package (v2.0.0)}

Check here (in particular the answer from Prof. Baddeley), the references therein, and the help page of st_sample for more details.