Poisson distribution sampling in 2D space with spatially varying intensity?

I am sampling points from a Poisson process in a 2D space (Ny, Nx) grid with a constant intensity (lam) using the code given below.

How can I sample values from Poisson process in a space where the intensity is not constant, but varying spatially/regionally?

import numpy as np
from scipy.spatial import cKDTree as kdtree
import matplotlib.pyplot as plt

Nx, Ny, n_cells_reject_criteria = 100, 100, 3
valid = False

while not valid:
    rate_lambda = 0.02
    #===========generate random samples from homogeneous poisson process===========
    mean_poisson = rate_lambda*Nx*Ny
    n_events_pp = np.random.poisson(lam=mean_poisson)
    x_pp = np.round(np.random.uniform(low=0, high=Nx-1, size=n_events_pp)) # generate n uniformly distributed points
    y_pp = np.round(np.random.uniform(low=0, high=Ny-1, size=n_events_pp)) # generate n uniformly distributed points
    coords_random_ji = ([np.int(j) for j in y_pp], [np.int(i) for i in x_pp])

    #===========test there are no adjacent cells===========
    valid = len(kdtree(coords_random_ji).query_pairs(n_cells_reject_criteria)) == 0

#===========plot resuls===========
#------- create an empty mesh
grid = np.zeros((Ny, Nx), dtype=np.bool)

#------- superimpose the results from rejection sampling
grid[coords_random_ji] = True

#------- create empty figure
fig = plt.figure(figsize=(5, 5)) # in inches
#------- plot
plt.imshow(grid)