sample points on polygon edge in R

Question

If I have a spatialpolygons object in R, how can I generate a set of n points that are on the edge of that polygon?

I originally thought I could just sample from the polygon vertices, but it looks like there are sometimes stretches where there are no vertices because the polygon edge is a straight line...

Answer 1

A simple solution would be to use st_segmentize() from package sf that adds points to straight lines, then sample along these finer points.

st_segmentize() has an argument dfMaxLength that defines the max distance to allow along a line. The smaller you will set, the more points you will have. It should at least be as small as the minimum distance between any two points.

library(sf)
library(tidyverse)

## original form
poly <- st_polygon(x=list(cbind(x=c(1,2,3,1),y=c(1,2,1,1))))

# segmentize, then convert to points
poly_points <- st_segmentize(poly, dfMaxLength = 0.1) %>% 
  st_coordinates() %>% 
  as.data.frame() %>% 
  select(X, Y) %>% 
  st_as_sf(coords = c("X", "Y"))

## plot: you can just use sample() now on your point dataset
plot(poly, reset = FALSE, main = "segmentize (black point), then sample 5 (red points)")
plot(poly_points, reset = FALSE, add = TRUE)
plot(poly_points[sample(1:nrow(poly_points), size = 5),], add = TRUE, col = 2, pch = 19)

To get the minimum distance between any two points (beware of zero):

 poly %>% 
  st_coordinates() %>% 
  as.data.frame() %>% 
  st_as_sf(coords = c("X", "Y")) %>% 
  st_distance() %>%  c() %>% 
  unique() %>% 
  sort

Answer 2

Assuming you want to draw points around the perimeter, I would split this into two parts:

P(Point p on Edge e) = P(point p | Edge e) P(Edge e)

with P(Edge e) proportional to its length. So first sample an edge, then sample a point on it.

Here's an example triangle:

  poly <- Polygon(list(x=c(1,2,3,1),y=c(1,2,1,1)))

We'll calculate the lengths of the sides:

  require(gsl) #for fast hypot function
  xy <- poly@coords
  dxy <- diff(xy)
  h <- hypot(dxy[,"x"], dxy[,"y"])

and draw a random side:

  e <- sample(nrow(dxy), 1, probs=h)

and then draw a point on that edge:

  u <- runif(1)
  p <- xy[e,] + u * dxy[e,]

Wrapping the whole thing in a function, we have:

  rPointOnPerimeter <- function(n, poly) {
    xy <- poly@coords
    dxy <- diff(xy)
    h <- hypot(dxy[,"x"], dxy[,"y"])

    e <- sample(nrow(dxy), n,replace=TRUE, prob=h)

    u <- runif(n)
    p <- xy[e,] + u * dxy[e,]

    p
  }

with a demo:

plot( rPointOnPerimeter(100,poly) )