Search code examples
rgisrastershapefileterra

Sampling random points in R with both raster and shapefile constraints

I am trying to sample random points according to 2 conditions:

  1. in non-NA cells of a raster
  2. inside the polygons of a shapefile

I tested the terra::spatSample() function, but I can't meet the second condition.

enter image description here

Here is a reproducible example:

library(terra)
set.seed(1)
f <- system.file("ex/elev.tif", package="terra")
f <- terra::rast(f)

v <- system.file("ex/lux.shp", package="terra")
v <- vect(v)
v <- subset(v, v$NAME_2 != "Mersch")

test <- terra::crop(f, v, mask = TRUE)

sp <- terra::spatSample(test, 1000, method="random", replace=FALSE, na.rm=TRUE, as.df=TRUE, as.points=TRUE, xy=TRUE, warn=TRUE, weights=NULL)

plot(test)
plot(v, add = TRUE)
plot(sp, add = TRUE)
Solution

  • One option is to convert f to a point SpatVector first, and then intersect() the resulting points with v. That way, only points that intersect with v will be returned. You can then sample() your points. Note that the use of set.seed() in this reprex means the sample() will not be 'random'.

    library(terra)
library(tidyterra)
library(ggplot2)

set.seed(1)
f <- system.file("ex/elev.tif", package="terra")
f <- rast(f)

# Convert f to points SpatVector and remove NAs
p <- as.points(f, na.rm = TRUE)

v <- system.file("ex/lux.shp", package="terra")
v <- vect(v)
v <- subset(v, v$NAME_2 != "Mersch")

# Return p if inside v
p1 <- intersect(p, v)

# Return difference between p and p1 (for illustrative purposes)
p2 <- erase(p, p1)

# Sample p1
sp <- sample(p1, 1000)

# Plot
ggplot() +
  geom_spatraster(data = f) +
  geom_spatvector(data = v, alpha = 0.6) +
  geom_spatvector(data = p2, colour = "#920000") +
  geom_spatvector(data = sp) +
  coord_sf(xlim = c(5.8, 5.9),
           ylim = c(49.5, 49.6)) +
  scale_fill_continuous(na.value = "transparent") +
  scale_x_continuous(breaks = seq(5.8, 5.9, length.out = 3)) +
  theme(panel.background = element_blank())

    result