1

I have a set of existing sf points and want to randomly sample 1 point from each of a 1km2 grid cell over a landscape.

I can create a grid of 1km2 polygons, but I can't work out in sf how to then sample points within these polygons.

E.g.

library("sf")
polys <- st_read(system.file("shape/nc.shp", package="sf"))
polys <- polys[1:3,]
polys <- st_geometry(polys) %>% st_cast("POLYGON")

pts <- st_sample(polys, 30)

plot(polys)
plot(pts, add = TRUE)

example points

Now I want to extract one element of pts for each polygon in polys, with the result looking like:

example outcome

Is there a function to do this within the sf framework?

1

Create a vector with the number of polygons:

spts <- numeric(length = length(polys))

Get a table of which points are in each polygon:

pwin <- t(st_within(pts, polys))
pwin
> Sparse geometry binary predicate list of length 3, where the predicate was `within'
   1: 2, 5, 9, 12, 14, 16, 21, 23, 25, 27
   2: 4, 7, 10, 11, 26, 28, 29
   3: 1, 3, 6, 8, 13, 15, 17, 18, 19, 20, ...

Randomly sample a single point in each polygon: (NB: the ifelse(is_empty(... call saves this from breaking where there are no points in a polygon, while second ifelse( save misleading failure where there is a single point in each polygon, see comment from @Spacedman)

library("purrr")

for(i in 1:length(polys)){

  spts[i] <- ifelse(is_empty(pwin[[i]]), 0, ifelse(length(pwin[[i]]) == 1, pwin[[i]][1], sample(pwin[[i]], 1)))
}

nspts <- pts[spts]

Viola

plot(polys)
plot(pts, add = TRUE)
plot(nspts, add = TRUE, col = "red", pch = 19)

enter image description here

  • 1
    This will fail badly if there is only one point in a polygon, because sample(3:5,1) samples from 3,4,5, and sample(4:5,1) samples from 4,5, but if you do sample(5,1) you get a sample from 1,2,3,4,5, and not 5 every time. See ?sample and use a safer sampling function such as resample defined in the help page. – Spacedman Jul 5 at 6:24
  • @Spacedman you are quite correct and for posterity I will add in my fix for this – Umaomamaomao Jul 5 at 6:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.