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