### NMST543 Spatial Statistics, exercise 9 ### ############################################## library(geoR) # Today we discuss the problem of spatial prediction (kriging, # named after Danie Gerhardus Krige). # The ultimate answer is that under the assumption of finite # second moments the optimal prediction (in the mean square sense) # is the conditional expectation E[ Z(x_0) | Z_n ]. # However, this is difficult to work with, except for # Gaussian random fields. Hence we look for linear predictions. # We study the same dataset as before: points(elevation,cex.min=1,cex.max=4) plot(elevation,lowess=TRUE) # Before the kriging itself we need to specify the covariance # structure of the random field. We estimate a parametric model # using MLE (see the previous exercise class). # Model with constant mean: ml0 <- likfit(elevation, ini=c(3000,2), cov.model="matern", kappa=1.5) ml0 # Model with linear trend: ml1 <- likfit(elevation, trend="1st", ini=c(1300,2), cov.model="matern", kappa=1.5) ml1 # Model quadratic trend: ml2 <- likfit(elevation, trend="2nd", ini=c(1300,2), cov.model="matern", kappa=1.5) ml2 ### Simple kriging ### ###################### # Allows non-constant mean value function but we have to know it # or at least estimate it. # Spatial grid for prediction: locs <- pred_grid(c(0,6.3), c(0,6.3), by=0.1) # Simple kriging with stationary model: KC <- krige.control(type="sk", obj.mod=ml0) # setting up the procedure sk <- krige.conv(elevation, krige=KC, loc=locs) # the estimation itself # Simple kriging with linear trend: KCt <- krige.control(type="sk", obj.mod=ml1, trend.d = "1st", trend.l="1st") skt <- krige.conv(elevation, krige=KCt, loc=locs) # Simple kriging with quadratic trend: KCt2 <- krige.control(type="sk", obj.mod=ml2, trend.d = "2nd", trend.l="2nd") skt2 <- krige.conv(elevation, krige=KCt2, loc=locs) # Range of values necessary for plotting the predictions: pred.lim <- range(c(sk$pred, skt$pred, skt$pred2)) pred.lim # The predictions do not leave the range of data values: summary(elevation) # For all plots today: dark color = high values; # light colors = low values. # Stationary model: windows() image(sk, col=gray(seq(1,0,l=51)), zlim=pred.lim) contour(sk, add=T, nlev=6) points(elevation, add=TRUE, cex.max=2) # Linear trend: windows() image(skt, col=gray(seq(1,0,l=51)), zlim=pred.lim) contour(skt, add=T, nlev=6) points(elevation, add=TRUE, cex.max=2) # Quadratic trend: windows() image(skt2, col=gray(seq(1,0,l=51)), zlim=pred.lim) contour(skt2, add=T, nlev=6) points(elevation, add=TRUE, cex.max=2) # Range of values necessary for plotting the standard # deviations of the predictions: sd.lim <- range(sqrt(c(sk$krige.var, skt$krige.var, skt2$krige.var))) # Stationary model: windows() image(sk, value=sqrt(sk$krige.var), col=gray(seq(1,0,l=51)), zlim=sd.lim) contour(sk, value=sqrt(sk$krige.var), levels=seq(10,27,by=2), add=T) points(elevation$coords, pch="+") # Linear trend: windows() image(skt, value=sqrt(skt$krige.var), col=gray(seq(1,0,l=51)), zlim=sd.lim) contour(skt, value=sqrt(skt$krige.var), levels=seq(10,27,by=2), add=T) points(elevation$coords, pch="+") # Quadratic trend: windows() image(skt2, value=sqrt(skt2$krige.var), col=gray(seq(1,0,l=51)), zlim=sd.lim) contour(skt2, value=sqrt(skt2$krige.var), levels=seq(10,27,by=2), add=T) points(elevation$coords, pch="+") # Higher variability (error) for models with linear or quadratic trend # is caused precisely by the need to estimate the parameters of the trend. # This introduces additional variability into the prediction. ### Ordinary kriging ### ######################## # We assume a constant mean value, the prediction is a linear combination # of the observed values. The sum of weights = 1; we allow negative # weights and weights > 1. # Simulated dataset (we need a new one, 'elevation' has some trend): summary(s100) points(s100) points(s100,pt.divide="quintile") plot(s100) # Estimate the covariance structure, using least # squares method now: bin1 <- variog(s100,uvec=seq(0,1,l=11)) plot(bin1) # Fit a model with zero nugget wls <- variofit(bin1,ini=c(1,0.5),fix.nugget=T) wls # Four locations for prediction: plot(s100$coords,xlim=c(0,1.2),ylim=c(0,1.2),xlab="X",ylab="Y") polygon(x=c(0,1,1,0),y=c(0,0,1,1),lty=2) loci <- matrix(c(0.2,0.6,0.2,1.1,0.2,0.3,1.0,1.1),ncol=2) text(loci, as.character(1:4), cex=1.3, col="red") # Prediction, ordinary kriging, parameters of variogram # are estimated in 'wls': kc4 <- krige.conv(s100,locations=loci,krige=krige.control(obj.mod=wls, type="ok")) # Predicted values: kc4$predict # Prediction errors: kc4$krige.var # Predictions in selected observation locations # (spatial prediction interpolates the data): kc3 <- krige.conv(s100,locations=s100$coords[1:3,],krige=krige.control(obj.m=wls, type="ok")) s100$data[1:3] kc3$predict # Spatial grid for prediction: pred.grid <- expand.grid(seq(0,1,l=51),seq(0,1,l=51)) # Prediction on the grid: kc <- krige.conv(s100,loc=pred.grid,krige=krige.control(obj.m=wls, type="ok")) image(kc,loc=pred.grid,col=gray(seq(1,0.1,l=30)),xlab="X",ylab="Y") points(s100,add=T) ### Universal kriging ### ######################### # Approach similar to ordinary kriging, we consider non-constant # trend given by a linear model. ### Individual task: ### ######################## # Try the universal kriging for the 'elevation' data. # Hint: check the help files for 'krige.conv', argument 'type.krige'.