How to get the same results using`rpois` and `rpoispp` from spatstat

I'm struggling using the functions rpois and rpoispp from the spatstat package, specifically when using lambda as a pixel image object. I want to sample points from certain pixels using lambda, but to get the same results as with rpois, I have to consider the number of cells lambda has. I guess this might be related to lambda vs intensity.

Why do I have to consider n.cells to use rpoispp in order to get the same results as with rpois? Is this the correct way to use rpoispp?

Here's an example:

library(spatstat)

set.seed(123)

# pixel image
Z <- as.im(function(x,y){10 * sqrt(x+y)}, unit.square())
n.cells <- dim(Z)[1]*dim(Z)[2]

# point pattern
pp <- rpoispp(lambda = Z) 
npoints(pp)
# [1] 8

# point pattern using n.cells
pp.cells <- rpoispp(lambda = Z*n.cells) 
npoints(pp.cells)
# [1] 160100

# counts
counts <- rpois(n=n.cells, lambda=Z[])
counts.im <- im(matrix(counts, nrow=dim(Z)[1], ncol=dim(Z)[2]))
sum(counts)
# [1] 159980

plot(Z)
plot(pp, add=T)

plot(pixellate(pp))

plot(Z)
plot(pp.cells, add=T)

plot(pixellate(pp.cells))

plot(counts.im)

The docs for rpoispp specifically say:

warning Note that lambda is the intensity, that is, the expected number of points per unit area. The total number of points in the simulated pattern will be random with expected value mu = lambda * a where a is the area of the window win.

If we examine Z, you will see it has an area of one square unit:

Z
#> real-valued pixel image
#> 128 x 128 pixel array (ny, nx)
#> enclosing rectangle: [0, 1] x [0, 1] units

This means that each cell has an area of 1/(128^2) or 1/16384 square units. Therefore, if we do:

rpoispp(lambda = Z) 

then for each cell, our lambda value will be the numerical value of the cell multiplied by 1/16384. Since this is at most 0.00086 for any individual cell, we will get no points drawn in almost every cell, with a small number of cells getting 1 point.

The expected value of the total number of points created by rpoispp(lambda = Z) is

sum(Z[])/(128^2)
#> [1] 9.75164

which is of course the same as

mean(Z[])
#> [1] 9.75164

If I understand correctly, you want to draw an average of n points for each cell, where n is equivalent to the cell's value. If this is the case, you need to correct for each cell's area.

In your example, this is straightforward, since each cell is 1/16384 square units in area, so we need only multiply the lambda of each cell by 16384 (i.e. n.cells) to get the same expected number of points we would with rpois(n.cells, Z[])

Therefore in this case, rpoispp(lambda = Z * n.cells) works, because your cells are in a 1-unit square, and they each have an area of 1/n.cells.

Perhaps a more natural way to do this without having to multiply by the number of cells is to rescale Z so that its cells each have an area of 1 unit:

set.seed(123)

Z <- as.im(function(x,y){10 * sqrt(x + y)}, unit.square())
Z <- rescale(Z, 1/128)
pp <- rpoispp(lambda = Z) 
npoints(pp)
#> [1] 159546

plot(Z)
plot(pp, add = TRUE, cex = 0.1)