# Algorithm for generating a Poisson process on a complicated 2d geometry

I am looking at some count data by geographic counties in California. As a starting point, a Poisson process came to mind--though there are other good choices like negative binomial, etc.

Given a $$\lambda$$ parameter for a county, I need to simulate a Poisson realization of events on that county geometry. Since the geometry of a county is rather irregular, I was wondering what a good algorithm might be to do something like this.

If I was using a simple rectangular county--for the sake of simplicity--then I was multiply the area of the rectangle $$A$$ times the intensity parameter $$\lambda$$, compute the number of events, and then generate that number of draws from a bivariate uniform distribution. That would give me the coordinates $$(length*x, height*y)$$ for the random event locations in a rectangle.

However if I have an irregular geometry, then the algorithm will generate points outside of the boundaries of a county.

I had a couple of ideas, but I was not sure if there were any unintended consequences of such approaches.

1. I could discretize the 2d county boundaries into a grid of coordinates. Then I would just have to pick the Poisson number of coordinates from this list randomly. Likely I would have to sample with replacement.

2. I could just sample points from a bounding box around the county boundaries. Then I could use the algorithm above to choose points in the bounding box. For the points that fall outside of the country boundaries, I would just count and discard those. Then I would try and again randomly assign the count of discarded points back to that bounding box. This time some of the new points will fall within the country boundaries and others would fall outside again. Repeat the process till all the points are assigned inside of the boundary.

These are just a couple of ideas. Can anyone tell me if one approach is better than the other. The question I am concerned with is how to ensure that I am still generating a poisson process and not introducing some unintentional dependencies through the generation process.

You can find the areas of the counties in California from https://en.wikipedia.org/wiki/List_of_counties_in_California. This will enable you to generate a count $$n_i$$ from the appropriate Poisson distribution (i.e., one with mean proportional to area) for each county $$i$$ with a trivial amount of effort. You can then apply a modified version of your algorithm 2 - generate uniformly distributed locations in the bounding box for each county, and rejecting those outside the county boundary - until the total number of locations inside the county boundary equals $$n_i$$.

This gets you Poisson counts uniformly distributed across each county with means proportional to the county areas, which seems to be what you want.

A very simple solution (in the same spirit as the answer by @jbowman) is to generate a Poisson point process in the bounding box and then restrict to the irregular region. This is valid since for a Poisson process on any region the restriction to any sub-region is also a Poisson process. So the steps are:

1. Find bounding box and calculate area
2. Calculate mean = lambda * area_of_box
3. Generate n from Poisson distribution with mean given above
4. Place the n points uniformly in box
5. Remove points outside the region of interest

The first four steps are computationally trivial while step 5 contains a point-in-polygon check which may be costly. However, you avoid having to make a triangulation.

Hmm, after some more research, it seems that the best way to go is triangulation. For a complex geometry, we can first triangular the geometry using any common meshing software like Gmsh or others. Once we triangulate the geometry, we can compute the area of that triangle and multiply that by the $$\lambda$$ parameter. Then we can generate the random Poisson points inside of each triangle.

Here is some code for computing the poisson process on a triangle, courtesy of this blog post.

#import libraries
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt

#Simulation window parameters -- points A,B,C of a triangle
xA=0;xB=1;xC=1; #x values of three points
yA=0;yB=0;yC=1; #y values of three points

#Point process parameters
lambda0=100; #intensity (ie mean density) of the Poisson process

#calculate sides of trinagle
a=np.sqrt((xA-xB)**2+(yA-yB)**2);
b=np.sqrt((xB-xC)**2+(yB-yC)**2);
c=np.sqrt((xC-xA)**2+(yC-yA)**2);
s=(a+b+c)/2; #calculate semi-perimeter

#Use Herron's forumula to calculate area -- or use polyarea
areaTotal=np.sqrt(s*(s-a)*(s-b)*(s-c)); #area of triangle

#Simulate a Poisson point process
numbPoints = scipy.stats.poisson(lambda0*areaTotal).rvs();#Poisson number of points
U = scipy.stats.uniform.rvs(0,1,((numbPoints,1)));#uniform random variables
V= scipy.stats.uniform.rvs(0,1,((numbPoints,1)));#uniform random variables

xx=np.sqrt(U)*xA+np.sqrt(U)*(1-V)*xB+np.sqrt(U)*V*xC;#x coordinates of points
yy=np.sqrt(U)*yA+np.sqrt(U)*(1-V)*yB+np.sqrt(U)*V*yC;#y coordinates of points

#Plotting
plt.scatter(xx,yy, edgecolor='b', facecolor='none', alpha=0.5 );
plt.xlabel("x"); plt.ylabel("y");