Search code examples
rrasterspatialterrar-stars

Sample each cell in a raster the number of points equal to the cell value

TLDR: I have a Poisson process raster and I want to turn it back into Log-Cox process raster

I have a raster which is an accumulation of counts in a cell. I need to sample as many points in each cell as the value of raster in each cell.

I can sample 1 observation from (hopefully) every cell like so:

library(terra)
#> terra 1.7.46
set.seed(42)
make_pois_raster <- function(x,y){rast(matrix(rpois(x*y, 2), x, y))}
r25 <- make_pois_raster(5,5)
plot(r25)
spatSample(r25, size=25, as.points=TRUE) |> plot(add=TRUE)
#> Warning: [is.lonlat] assuming lon/lat crs

What I would really like to do is something like this (does not currently work)

spatSample(r25, size=values(r25, mat=FALSE), as.points=TRUE)

spatSample can sample multiple points per geometry if sampling from a vector, but for raster only a fixed number of points can be sampled.

Created on 2023-09-13 with reprex v2.0.2

Solution

  • Could you perhaps convert the raster to polygons and sample from those instead using the size parameter?

    # convert raster to polygons
v <- terra::as.polygons(r25, aggregate = FALSE)

# remove polygons having value 0
v <- terra::subset(v, v$lyr.1 > 0)

plot(v)

spatSample(x, size = unlist(terra::values(v))) |> plot(add = TRUE)

    enter image description here

    Note the use of aggregate = FALSE when converting to polygons, and also that spatSample() has a problem with the size vector including zeros, so these are removed.