#### Bayesian kriging and spatial prediction:Surface elevation

The data for this example are taken from Diggle and Riberio (2000) (and came originally from Davis (1973)). The data file contains the variables height, x and y, giving surface elevation at each of 52 locations (x, y) within a 310-foot square. The unit of distance is 50 feet, whereas one unit of height represents 10 feet of elevation. A Gaussian kriging model can be fitted to these data in OpenBUGS using either the spatial.exp or spatial.disc distributions. The data file also contains a set of locations x.pred and y.pred representing a 15 x 15 grid of points at which we wish to predict surface elevation. Predictions can be obtained using either the spatial.pred or spatial.unipred predictive distributions in OpenBUGS

Model
model {
# Spatially structured multivariate normal likelihood # exponential correlation function
height[1:N] ~ spatial.exp(mu[], x[], y[], tau, phi, kappa)
# disc correlation function
#   height[1:N] ~ spatial.disc(mu[], x[], y[], tau, alpha)

for(i in 1:N) {
mu[i] <- beta
}
# Priorsbeta ~ dflat()
tau ~ dgamma(0.001, 0.001)
sigma2 <- 1/tau

# priors for spatial.exp parameters
# prior range for correlation at min distance (0.2 x 50 ft) is 0.02 to 0.99
phi ~ dunif(0.05, 20)
# prior range for correlation at max distance (8.3 x 50 ft) is 0 to 0.66
kappa ~ dunif(0.05,1.95)

# priors for spatial.disc parameter
# prior range for correlation at min distance (0.2 x 50 ft) is 0.07 to 0.96
#   alpha ~ dunif(0.25, 48)
# prior range for correlation at max distance (8.3 x 50 ft) is 0 to 0.63
# Spatial prediction
# Single site prediction
for(j in 1:M) {
height.pred[j] ~ spatial.unipred(beta, x.pred[j], y.pred[j], height[])
}

# Only use joint prediction for small subset of points, due to length of time it takes to run
for(j in 1:10) { mu.pred[j] <- beta }
height.pred.multi[1:10] ~ spatial.pred(mu.pred[], x.pred[1:10], y.pred[1:10], height[])
}
}

##### Data
``` list(N=52, M= 225, x = c(0.3, 1.4, 2.4, 3.6, 5.7, 1.6, 2.9, 3.4, 3.4, 4.8, 5.3, 6.2,    0.2, 0.9, 2.3, 2.5, 3, 3.5, 4.1, 4.9, 6.3, 0.9, 1.7, 2.4, 3.7, 4.5, 5.2,   6.3, 0.3, 2, 3.8, 6.3, 0.6, 1.5, 2.1, 2.1, 3.1, 4.5, 5.5, 5.7, 6.2, 0.4,   1.4, 1.4, 2.1, 2.3, 3.1, 4.1, 5.4, 6, 5.7, 3.6),    y = c(6.1, 6.2, 6.1, 6.2, 6.2, 5.2, 5.1, 5.3, 5.7, 5.6, 5, 5.2, 4.3,    4.2, 4.8, 4.5, 4.5, 4.5, 4.6, 4.2, 4.3, 3.2, 3.8, 3.8, 3.5, 3.2, 3.2,    3.4, 2.4, 2.7, 2.3, 2.2, 1.7, 1.8, 1.8, 1.1, 1.1, 1.8, 1.7, 1, 1, 0.5,    0.6, 0.1, 0.7, 0.3, 0, 0.8, 0.4, 0.1, 3, 6),    height = c(870, 793, 755, 690, 800, 800, 730, 728,   710, 780, 804, 855, 830, 813, 762, 765, 740, 765, 760, 790, 820, 855,    812, 773, 812, 827, 805, 840, 890, 820, 873, 875, 873, 865, 841, 862,    908, 855, 850, 882, 910, 940, 915, 890, 880, 870, 880, 960, 890, 860,    830, 705),   x.pred=c(0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0.21,    0.21, 0.21, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63,    0.63, 0.63, 0.63, 0.63, 0.63, 1.05, 1.05, 1.05, 1.05, 1.05, 1.05, 1.05,    1.05, 1.05, 1.05, 1.05, 1.05, 1.05, 1.05, 1.05, 1.47, 1.47, 1.47, 1.47,    1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.89,    1.89, 1.89, 1.89, 1.89, 1.89, 1.89, 1.89, 1.89, 1.89, 1.89, 1.89, 1.89,    1.89, 1.89, 2.31, 2.31, 2.31, 2.31, 2.31, 2.31, 2.31, 2.31, 2.31, 2.31,    2.31, 2.31, 2.31, 2.31, 2.31, 2.73, 2.73, 2.73, 2.73, 2.73, 2.73, 2.73,    2.73, 2.73, 2.73, 2.73, 2.73, 2.73, 2.73, 2.73, 3.15, 3.15, 3.15, 3.15,    3.15, 3.15, 3.15, 3.15, 3.15, 3.15, 3.15, 3.15, 3.15, 3.15, 3.15, 3.57,    3.57, 3.57, 3.57, 3.57, 3.57, 3.57, 3.57, 3.57, 3.57, 3.57, 3.57, 3.57,    3.57, 3.57, 3.99, 3.99, 3.99, 3.99, 3.99, 3.99, 3.99, 3.99, 3.99, 3.99,    3.99, 3.99, 3.99, 3.99, 3.99, 4.41, 4.41, 4.41, 4.41, 4.41, 4.41, 4.41,    4.41, 4.41, 4.41, 4.41, 4.41, 4.41, 4.41, 4.41, 4.83, 4.83, 4.83, 4.83,    4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 4.83, 5.25,    5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25,    5.25, 5.25, 5.67, 5.67, 5.67, 5.67, 5.67, 5.67, 5.67, 5.67, 5.67, 5.67,    5.67, 5.67, 5.67, 5.67, 5.67, 6.09, 6.09, 6.09, 6.09, 6.09, 6.09, 6.09,    6.09, 6.09, 6.09, 6.09, 6.09, 6.09, 6.09, 6.09), y.pred=c(6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31, 1.89, 1.47, 1.05,    0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31,    1.89, 1.47, 1.05, 0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57,    3.15, 2.73, 2.31, 1.89, 1.47, 1.05, 0.63, 0.21, 6.09, 5.67, 5.25, 4.83,    4.41, 3.99, 3.57, 3.15, 2.73, 2.31, 1.89, 1.47, 1.05, 0.63, 0.21, 6.09,    5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31, 1.89, 1.47, 1.05,    0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31,    1.89, 1.47, 1.05, 0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57,    3.15, 2.73, 2.31, 1.89, 1.47, 1.05, 0.63, 0.21, 6.09, 5.67, 5.25, 4.83,    4.41, 3.99, 3.57, 3.15, 2.73, 2.31, 1.89, 1.47, 1.05, 0.63, 0.21, 6.09,    5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31, 1.89, 1.47, 1.05,    0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31,    1.89, 1.47, 1.05, 0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57,    3.15, 2.73, 2.31, 1.89, 1.47, 1.05, 0.63, 0.21, 6.09, 5.67, 5.25, 4.83,    4.41, 3.99, 3.57, 3.15, 2.73, 2.31, 1.89, 1.47, 1.05, 0.63, 0.21, 6.09,    5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31, 1.89, 1.47, 1.05,    0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57, 3.15, 2.73, 2.31,    1.89, 1.47, 1.05, 0.63, 0.21, 6.09, 5.67, 5.25, 4.83, 4.41, 3.99, 3.57,    3.15, 2.73, 2.31, 1.89, 1.47, 1.05, 0.63, 0.21)) ```

Initial values for exponential model

##### Inits for chain 1
``` list(tau=0.001, phi=0.4, beta=820, kappa=1,height.pred=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),height.pred.multi=c(0,0,0,0,0,0,0,0,0,0)) ```

##### Inits for chain 2
``` list(tau=0.1, phi=0.1, beta=1000, kappa=0.1,height.pred=c(100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100, 100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100),height.pred.multi=c(100,100,100,100,100,100,100,100,100,100)) ```

Initial values for disc model

##### Inits for chain 1
``` list(tau=0.001, phi = 0.1, beta=820, alpha=5,height.pred=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),height.pred.multi=c(0,0,0,0,0,0,0,0,0,0)) ```

##### Inits for chain 2
``` list(tau=0.1, phi=0.1, beta=1000, alpha = 1,height.pred=c(100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100, 100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100),height.pred.multi=c(100,100,100,100,100,100,100,100,100,100)) ```

The GeoBUGS map tool can be used to produce an approximate map of the posterior mean predicted surface elevation. It is not possible to map point locations using GeoBUGS. However, a map file called elevation is already loaded in the GeoBUGS Map Tool --- this is a 15 x 15 grid with the centroids of the grid cell corresponding to the locations of each of the prediction points. You can use this to produce a map of the posterior mean (or other posterior summaries) of the vector of predicted values height.pred.

Results for exponential model

Results for disc model