Generating random points in polygons to simulate population distribution

I have a spatialpolygondataframe that represents census tracts. I want to simulate different scenarios of population density or distribution to test a model. To do this, I want to generate random points within the census tracts to represent the population. However, I want to skew the location of the population density or distribution. For example, in scenario 1 the location of the highest population density is west and the location of high population density is east in scenario 2.

I found an answer to a similar question here however the point distribution is skewed based on population size. Thus, I can't use the solution for my simulation. I tried the mock code below based on the solution

library(raster)
library(sf)
library(dplyr)

nsplit = function(X,n){
p = X/sum(X)
diff(round(n*cumsum(c(0,p))))
}

st_transform(32617)

sample_points <- st_sample(nc, size = 1000,
type = 'random', exact = TRUE) %>%
# optional: Give the points an ID
st_sf('ID' = seq(length(.)), 'geometry' = .) %>%
# optional: Get underlying polygon attributes
st_intersection(., nc)

plot(nc, reset = FALSE, axes = TRUE)
plot(sample_points['NAME'], pch = 19, cex = 0.5, add = TRUE)

And this is the point distribution pattern I got Using the NC polygons in the figure above, I would like to have a high number of points in the west for my first simulation scenario and east of NC for the second.

• What you want is a non-uniform point pattern. You have to decide exactly how non-uniform you want it, by defining an intensity function, which is a function of x and y (or lat and long). All you've said is that you want the intensity to increase or decrease with X coordinate, but not by how much, or if it also needs to relate to population, or if you want it linear or not. All these things need to be defined before we can give a complete solution. Nov 4 '21 at 18:57
• @Spacedman this field is completely new to me hence the limited amount of information on the issues you raised. I don't want the intensity function to relate to the population. The purpose of the simulate is to see if the result from the model which seems to show actual population distribution is real or artefactual. For now, a linear intensity function will suffice and any magnitude for the increase or decrease function will suffice (I can later change these for further simulation analysis if the need be) Nov 4 '21 at 19:14
• If you want to test a statistical model you generally want to simulate from that model's distribution. Is that what you have? Nov 4 '21 at 19:57
• You do not want a random sample but, rather a non stationary point process. If you think about it, density is not random, its clustered. As @Spacedman stated you can do this via an intensity function but, you really need to decide on an underlying distributional form. However, it would be ideal if you had a model, based on observed data, to start with. Else wise, the results are rather arbitrary in their meaning. I may help if stated something resembling a hypothesis, regarding the spatial process, that you want to test. Nov 4 '21 at 22:29
• I am currently analysing disease distribution (Bayesian spatial modeling) at a point level (household) accounting for some underlying spatial and individual covariates. Now I want to create fake data (with the ability to change the population distribution and disease distribution) to test the model. I will then run the model with the stimulated data, to see if the model is just mimicking the underlying population distribution or the actual disease pattern. It is also the reason for creating simulated data of different scenarios. I apologise if my submissions are not helpful - a new field to me Nov 5 '21 at 6:14

Here's a way to do something via rejection sampling. There's probably a function to do this somewhere, or a better way to do it in the linear case, but writing it from scratch isn't that hard, plus it makes it easier to replace the linear population density gradient with any other function if you want.

This function simulates non-uniform random points in a bounding box defined by xmin,xmax and ymin, ymax. The points are uniform in Y coordinate, but non-uniform in X coordinate. The density is linear across Y and proportional to z0 at xmin and z1 at xmax.

ewpop <- function(N, xmin, xmax, ymin, ymax, z0, z1){

x = runif(N, xmin, xmax)

## linear density function
zf = z0 + (z1-z0)*(x-xmin)/(xmax-xmin)

## rejection probability
p = runif(N, 0, max(c(z0,z1)))

## accept these:
x = x[p<zf]

## make some Y coords for those:
y = runif(length(x), ymin, ymax)

cbind(x,y)
}

Example:

Generate some points in (0,3) (0,1). The histogram of the X coordinate shows the non-uniform distribution.

> p1 = ewpop(10000, 0,3,0,1,10,20)
> plot(p1, asp=1, pch=".")
> hist(p1[,1]) For a more extreme inhomogeneity, set z0 to 0:

> p1 = ewpop(10000, 0,3,0,1,0,1)
> plot(p1,asp=1,pch=".")
> hist(p1[,1]) Note that because this is rejection sampling you don't get back the number of points you asked for. In the worst case (which is when z0==0) you'll get approximately half.

> nrow(p1)
 5076

In which case just call it again and rbind everything together - each call is independent.

Anyway, to get points in your irregular study area you are going to have to intersect your shape with the points in the bounding box from this function and iterate until you have enough points. Or generate way too many with this function, then do the intersection, and if you have more than you need the only take the first N that you need - again they are all independent.

• Thank you, it works prefectly! Nov 5 '21 at 10:10