Doubly stochastic Bayesian hierarchical models are introduced to account for uncertainty and spatial variation in the underlying intensity measure for point process models. Inhomogeneous gamma process random fields and, more generally, Markov random fields with infinitely divisible distributions are used to construct positively autocorrelated intensity measures for spatial Poisson point processes; these in turn are used to model the number and location of individual events. A data augmentation scheme and Markov chain Monte Carlo numerical methods are employed to generate samples from Bayesian posterior and predictive distributions. The methods are developed in both continuous and discrete settings, and are applied to a problem in forest ecology. Bayesian mixture model; Bioabundance; Cox process; Data augmentation; Levy process; Markov chain Monte Carlo; Simulation.