The Gaussian geostatistical model has been widely used in spatial data modeling. In spite of its popularity, this model suffers from a severe implementation problem for Bayesian inference, for which a covariance matrix needs to be inverted at each iteration. This is infeasible when the number of observations is large. In this paper, we propose a predictionoriented Bayesian site selection (BSS) approach to tackle this difficulty. By dividing the observations into two sets, response variables and explanatory variables, the BSS approach forms a regression model which relates the observations through a conditional likelihood derived from the original Gaussian geostatistical model, and then reduces the dimension of the data using a stochastic variable selection procedure. Our numerical results indicate that the BSS approach can produce very good parameter estimates and prediction for large spatial data, while significantly reducing the computational time required by conventional Bayesian approaches.