  We propose a three-step divide-and-conquer strategy within the Bayesian
paradigm that delivers massive scalability for any spatial process model. We
partition the data into a large number of subsets, apply a readily available
Bayesian spatial process model on every subset, in parallel, and optimally
combine the posterior distributions estimated across all the subsets into a
pseudo-posterior distribution that conditions on the entire data. The combined
pseudo posterior distribution replaces the full data posterior distribution for
predicting the responses at arbitrary locations and for inference on the model
parameters and spatial surface. Based on distributed Bayesian inference, our
approach is called "Distributed Kriging" (DISK) and offers significant
advantages in massive data applications where the full data are stored across
multiple machines. We show theoretically that the Bayes $L_2$-risk of the DISK
posterior distribution achieves the near optimal convergence rate in estimating
the true spatial surface with various types of covariance functions, and
provide upper bounds for the number of subsets as a function of the full sample
size. The model-free feature of DISK is demonstrated by scaling posterior
computations in spatial process models with a stationary full-rank and a
nonstationary low-rank Gaussian process (GP) prior. A variety of simulations
and a geostatistical analysis of the Pacific Ocean sea surface temperature data
validate our theoretical results.
