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.

*

*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.


*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.

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");

