Overview

This document demonstrates how to generate random sampling points within a polygon while enforcing a minimum distance between points using a Simple Sequential Inhibition (SSI) process.

Libraries

library(sf)
library(spatstat.geom)
library(spatstat.random)
library(dplyr)
library(ggplot2)
library(maptiles)
library(ggspatial)
library(units)

We load:

• sf for modern spatial data handling
• spatstat modules for generating spatial point patterns
• ggplot2 for plotting
• maptiles and ggspatial for modern basemaps and map annotations
• units to safely handle distance units (e.g. metres)

First, create a bounding polygon using longitude and latitude coordinates in the WGS84 coordinate system. While this format is convenient for defining locations, it is not suitable for measuring distances because the units are in degrees rather than metres. To address this, the polygon is reprojected into the British National Grid coordinate system (EPSG:27700), which uses metres as its unit. This transformation is essential for ensuring that the minimum distance constraint applied later is accurate and interpretable. You can of course use our own shapefile if you don’t want to manually create a bounding polygon.

spatstat uses its own internal representation of space called an “observation window” (owin). We convert the sf polygon into owin format so it can be used as the boundary within which points will be generated. Without this conversion, the simulation functions in spatstat would not recognise the study area.

x <- c(-3.541238, -3.547341, -3.561397, -3.572132, -3.572495, -3.551687,
       -3.566528, -3.599789, -3.592162, -3.604156, -3.629145, -3.618190,
       -3.585571, -3.572779, -3.525398, -3.504483, -3.525454, -3.541238)

y <- c(51.68106, 51.69050, 51.7003, 51.69617, 51.68800, 51.67405,
       51.64949, 51.65832, 51.67560, 51.68177, 51.71373, 51.72463,
       51.73346, 51.73788, 51.71745, 51.69461, 51.67548, 51.68106)

# Create sf polygon in WGS84
poly_wgs <- st_polygon(list(cbind(x, y))) |> 
  st_sfc(crs = 4326)

# Transform to British National Grid (meters)
poly <- st_transform(poly_wgs, 27700)

win <- as.owin(poly)

Points are generated using the Simple Sequential Inhibition (SSI) process, which places points randomly but rejects any candidate point that falls within a specified distance of an existing point. In this case, the inhibition distance is set to 250 metres. The result is initially returned as a ppp object (a spatstat point pattern), which is then converted into an sf object for compatibility with the rest of the workflow.

samp1_ppp <- rSSI(r = 250, n = 10, win = win)  # 250 meters

# Convert to sf
samp1 <- st_as_sf(as.data.frame(samp1_ppp), coords = c("x", "y"), crs = 27700)

We then generate a full matrix of pairwise distances between all generated points. Inspecting this matrix provides a direct way to verify that the minimum distance constraint has been respected across all points.

dist_matrix <- st_distance(samp1)
dist_matrix

## Units: [m]
##           1        2        3         4         5        6        7        8
## 1     0.000 2135.996 4802.092 5196.1583 5886.0995 2902.544 3236.065 5579.278
## 2  2135.996    0.000 2686.834 3167.0465 3837.2100 3020.327 1567.009 4052.468
## 3  4802.092 2686.834    0.000 1781.4970 2072.7938 5059.138 2016.413 2797.764
## 4  5196.158 3167.047 1781.497    0.0000  695.3981 4397.598 3343.561 4577.688
## 5  5886.100 3837.210 2072.794  695.3981    0.0000 5030.483 3868.416 4783.126
## 6  2902.544 3020.327 5059.138 4397.5979 5030.4835    0.000 4579.629 7036.075
## 7  3236.065 1567.009 2016.413 3343.5610 3868.4157 4579.629    0.000 2494.419
## 8  5579.278 4052.468 2797.764 4577.6882 4783.1257 7036.075 2494.419    0.000
## 9  3258.136 2544.627 4093.057 3264.5085 3884.8693 1152.752 3989.559 6329.889
## 10 3566.796 2330.221 3392.810 2467.5012 3090.4724 1940.410 3587.215 5791.073
##            9        10
## 1  3258.1357 3566.7957
## 2  2544.6266 2330.2210
## 3  4093.0572 3392.8101
## 4  3264.5085 2467.5012
## 5  3884.8693 3090.4724
## 6  1152.7523 1940.4100
## 7  3989.5591 3587.2155
## 8  6329.8888 5791.0727
## 9     0.0000  797.1189
## 10  797.1189    0.0000

To summarise the spacing between points, the minimum non-zero distance for each point is extracted from the distance matrix. This represents the distance from each point to its nearest neighbour. The resulting values are added back to the spatial dataset as a new attribute (min_dist_m), providing a convenient way to assess how close each point is to its closest neighbour.

nearest <- apply(dist_matrix, 1, function(v) min(v[v > 0]))

samp1 <- samp1 |> 
  mutate(min_dist_m = as.numeric(set_units(nearest, "m")))

Web-based map tiles are typically served in geographic coordinates (WGS84), so the polygon is temporarily transformed back into that coordinate system for tile retrieval. The maptiles package is used to download a clean basemap from the CartoDB Positron provider, which is well suited for analytical visualisations due to its minimal styling.

The final map combines the basemap tiles with the study polygon and sampled points. Additional cartographic elements such as a scale bar and north arrow are added using ggspatial.

# Transform back to WGS84 for tile download
poly_wgs <- st_transform(poly, 4326)

tiles <- get_tiles(poly_wgs, provider = "CartoDB.Positron", zoom = 12)

# Plot 
ggplot() +
  layer_spatial(tiles) +
  geom_sf(data = poly, fill = "red", alpha = 0.3) +
  geom_sf(data = samp1, color = "blue", size = 2) +
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "tl", which_north = "true") +
  theme_minimal()

This function encapsulates the full workflow into a reusable tool. It takes a polygon, a minimum distance, and a desired number of points as inputs, then generates a spatial point dataset that satisfies the distance constraint. It also computes and attaches the nearest neighbour distance for each point.

# Function to generate random points within a polygon buffered by a minimum distance
# Attaches a dataframe of distances-between-points
#
# d = inherited distance (see rSSI documentation)
# n = number of points
# p = polygon
#
# E.g.
# spat <- Rpoints(d = 0.05, n = 50, p = poly1)
Rpoints <- function(d = 250, n = 10, p){
  
  stopifnot(inherits(p, "sf"))
  
  # Ensure projected CRS (meters)
  if (st_is_longlat(p)) {
    stop("Polygon must be projected (e.g., EPSG:27700)")
  }
  
  win <- as.owin(p)
  
  pts_ppp <- rSSI(r = d, n = n, win = win)
  
  pts <- st_as_sf(as.data.frame(pts_ppp), coords = c("x", "y"), crs = st_crs(p))
  
  dist_matrix <- st_distance(pts)
  
  nearest <- apply(dist_matrix, 1, function(v) min(v[v > 0]))
  
  pts |> 
    mutate(min_dist_m = as.numeric(set_units(nearest, "m")))
}

# Example usage
spat <- Rpoints(d = 200, n = 30, p = poly)
head(spat)