The paper presents an extension of geostatistical methods to situations where the data exhibit non-Gaussian variation. Traditional geostatistical methods, such as kriging, assume that the data are generated by a Gaussian spatial stochastic process. However, in some applications, such as assessing residual contamination from nuclear weapons testing on a South Pacific island and analyzing spatial variation in the risk of campylobacter infections, Gaussian assumptions are inappropriate. In the first example, the data are modeled as Poisson counts, conditional on an unobserved spatially varying intensity of radioactivity. In the second example, the data are modeled as binomial counts, conditional on an unobserved relative risk surface. The authors propose a Bayesian framework, implemented via Markov chain Monte Carlo (MCMC) methods, to solve the prediction problem for non-linear functionals of the spatial process. The theoretical framework is based on generalized linear mixed models, where the spatial process is treated as a Gaussian process and the data are modeled as a generalized linear model with the spatial process as an offset term. The authors demonstrate that the Bayesian approach allows for proper accounting of uncertainty in the estimation of model parameters. The paper also discusses the extension of the variogram to the generalized linear prediction model and presents applications to a simulated case-study and the two motivating examples. The results show that the Bayesian approach provides more accurate predictions and better accounts for uncertainty compared to traditional geostatistical methods.The paper presents an extension of geostatistical methods to situations where the data exhibit non-Gaussian variation. Traditional geostatistical methods, such as kriging, assume that the data are generated by a Gaussian spatial stochastic process. However, in some applications, such as assessing residual contamination from nuclear weapons testing on a South Pacific island and analyzing spatial variation in the risk of campylobacter infections, Gaussian assumptions are inappropriate. In the first example, the data are modeled as Poisson counts, conditional on an unobserved spatially varying intensity of radioactivity. In the second example, the data are modeled as binomial counts, conditional on an unobserved relative risk surface. The authors propose a Bayesian framework, implemented via Markov chain Monte Carlo (MCMC) methods, to solve the prediction problem for non-linear functionals of the spatial process. The theoretical framework is based on generalized linear mixed models, where the spatial process is treated as a Gaussian process and the data are modeled as a generalized linear model with the spatial process as an offset term. The authors demonstrate that the Bayesian approach allows for proper accounting of uncertainty in the estimation of model parameters. The paper also discusses the extension of the variogram to the generalized linear prediction model and presents applications to a simulated case-study and the two motivating examples. The results show that the Bayesian approach provides more accurate predictions and better accounts for uncertainty compared to traditional geostatistical methods.