c173c273_lec11_w11[1]

c173c273_lec11_w11[1] - University of California, Los...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: University of California, Los Angeles Department of Statistics Statistics C173/C273 Instructor: Nicolas Christou Ordinary kriging using geoR and gstat In this document we will discuss kriging using the R packages geoR and gstat. We will use the numerical example from last lecture. Here it is: A simple example: Consider the following data si s1 s2 s3 s4 s5 s6 s7 s0 x 61 63 64 68 71 73 75 65 y 139 140 129 128 140 141 128 137 z(si ) 477 696 227 646 606 791 783 ??? Our goal is to estimate the unknown value at location s0 . Here is the x - y plot: q s6 128 130 132 134 136 138 140 q q s2 q s5 s1 q s0 y coordinate q s3 q s4 70 72 s7 74 q 62 64 66 68 x coordinate 1 For these data, let's assume that we use the exponential semivariogram model with parameters c0 = 0, c1 = 10, = 3.33. (h) = c0 + c1 (1 - e- ) = 10(1 - e- 3.33 ). which is equivalent to the covariance function C(h) = c0 + c1 , h c1 e - , h=0 h>0 C(h) = 10, h 10e- 3.33 , h=0 h>0 h h The predicted value at location s0 is equal to: n z (s0 ) = ^ i=1 wi z(si ) = 0.174(477) + + 0.086(783) = 592.59. And the variance: n 2 e = i=1 wi (si - s0 ) + = 0.174(7.384) + + 0.086(9.823) + 0.906 = 8.96. Kriging using geoR: We will use now the geoR package to find the same result. First we read our data as a geodata object: > a <- read.table("http://www.stat.ucla.edu/~nchristo/statistics_c173_c273/ kriging_11.txt", header=TRUE) > b <- as.geodata(a) To predict the unknown value at locaton (x = 65, y = 137) we use the following: > prediction <- ksline(b, cov.model="exp", cov.pars=c(10,3.33), nugget=0, locations=c(65,137)) where, b The geodata cov.model The model we are using cov.pars The parameters of the model (partial sill and range) nugget The value of the nugget effect locations The coordinates (x, y) of the points to be predicted The object "prediction" contains among other things the predicted value at location x = 65, y = 137 and its variance. We can obtain them as follows: > prediction$predict [1] 592.7587 > prediction$krige.var [1] 8.960294 2 Suppose now we want to predict the value at many locations. The following commands will produce a grid whose points will be estimated using kriging: > > > > > x <- seq(61, 75, by=1) y <- seq(128,141, by=1) xv <- rep(x,14) yv <- rep(y, each=15) in_mat <- as.matrix(cbind(xv,yv)) Of course the matrix can be constructed also as follows: > > > > > > > > > > x.range <- as.integer(range(a[,1])) y.range <- as.integer(range(a[,2])) x=seq(from=x.range[1], to=x.range[2], by=1) y=seq(from=y.range[1], to=y.range[2], by=1) length(x) length(y) xv <- rep(x,length(y)) yv <- rep(y, each=length(x)) in_mat <- as.matrix(cbind(xv,yv)) plot(in_mat) This grid consists of 15 14 = 210 points stored in the matrix in mat. The command that predicts the value at each one of these points is the following: > q <- ksline(b, cov.model="exp",cov.pars=c(10,3.33), nugget=0, locations=in_mat) ksline: kriging location: 1 out of 210 ksline: kriging location: 101 out of 210 ksline: kriging location: 201 out of 210 ksline: kriging location: 210 out of 210 Kriging performed using global neighbourhood We can access the predicted values and their variances using q$predict and q$krige.var. Here are the first 5 predicted values with their variances: > cbind(q$predict[1:5], q$krige.var[1:5]) [,1] [,2] [1,] 458.4491 9.245493 [2,] 413.2103 7.850838 [3,] 362.4674 5.927999 [4,] 338.9828 4.516906 [5,] 393.3933 5.280417 3 To construct the raster map we type: image(q, val=q$predict) Or simply: image(q) Here is the plot: 140 Y Coord 128 130 132 134 136 138 65 X Coord 70 75 And here is the plot with the data points: points(a) q 140 q q q Y Coord 130 132 134 136 138 q 128 q q 65 X Coord 70 75 4 We can construct a raster map of the variances: > image(q, val=q$krige.var) Here is the plot: 140 Y Coord 128 130 132 134 136 138 65 X Coord 70 75 And also we can construct a raster map of the standard errors: > image(q, val=sqrt(q$krige.var)) Here is the plot: 140 Y Coord 128 130 132 134 136 138 65 X Coord 70 75 5 The following command will construct a perspective plot of the predicted values: > persp(x,y,matrix(q$predict,15,14), xlab="x coordinate", ylab="y coordinate", zlab="Predicted values of z", main="Perspective plot of the predicted values") Perspective plot of the predicted values yco ord ina te And here is the perspective plot of the standard errors: > persp(x,y,matrix(sqrt(q$krige.var),15,14), xlab="x coordinate", ylab="ycoordinate", zlab="Predicted values of z", main="Perspective plot of the standard errors") Perspective plot of the standard errors yco ord ina te lue Predictedva s ofz Predicted va lues of z x coordinate x coordinate 6 Commnets: The argument locations for the function ksline can be a matrix or a data frame. So far we use a matrix (in mat). We can also use a data frame as follows: > x.range <- as.integer(range(a[,1])) > y.range <- as.integer(range(a[,2])) > grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=1), y=seq(from=y.range[1], to=y.range[2], by=1)) > q <- ksline(b, cov.model="exp",cov.pars=c(10,3.33), nugget=0, locations=grd) Another function in geoR that performs kriging is the krige.conv function. It can be used as follows: > kc <- krige.conv(b, loc=in_mat, krige=krige.control(cov.pars=c(10, 3.33), nugget=0)) krige.conv: model with constant mean krige.conv: Kriging performed using global neighbourhood Or using a data frame argument for locations: > kc <- krige.conv(b, loc=grd, krige=krige.control(cov.pars=c(10, 3.33), nugget=0)) krige.conv: model with constant mean krige.conv: Kriging performed using global neighbourhood In case you have a variofit output you can use it as an input of the argument krige as follows (this is only an example): var1 <- variog(b, max.dist=1000) fit1 <- variofit(var1, cov.model="exp", ini.cov.pars=c(1000, 100), fix.nugget=FALSE, nugget=250) qq <- krige.conv(b, locations=grd, krige=krige.control(obj.model=fit1)) 7 Kriging using gstat: We will use now the gstat package to find the same result. First we read our data and create the grid for prediction as follows: > a <- read.table("http://www.stat.ucla.edu/~nchristo/statistics_c173_c273/ kriging_11.txt", header=TRUE) > x.range <- as.integer(range(a[,1])) > y.range <- as.integer(range(a[,2])) > grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=1), y=seq(from=y.range[1], to=y.range[2], by=1)) We now define the model. Normally the model must be estimated from the sample variogram, but for this simple example we assume that it is given as below: > library(gstat) > m <- vgm(10, "Exp", 3.33, 0) There are two ways to perform ordinary kriging with gstat. The data and the grid are used as data frames, with the extra argument locations as shown below: > q1 <- krige(id="z", formula=z~1, data=a, newdata=grd, model = m, locations=~x+y) The other way is to convert the data and the grid as spatial data points data frame: > coordinates(a) <- ~x+y > coordinates(grd) <- ~x+y > q2 <- krige(id="z", formula=z~1, a, newdata=grd, model = m) Important note: If we use the second way the argument data= is not allowed. We simply use the name of the data, here just a. Also, q1 is a data frame, while q2 is spatial data points data frame. Using q1 we can create a 3D plot with the libraries scatterploted and rgl as follows: > library(scatterplot3d) > library(rgl) > scatterplot3d(q1$x, q1$y, q1$z.pred, xlab="x", ylab="y", zlab="Predicted values") > plot3d(q1$x, q1$y, q1$z.pred, size=3) 8 Here are the plots: (a). Using the scatterplot3d command. q q qq q q q q q q q q q qq q q q q q q q Predicted values q q q q q q q q q q q q qq qq q q q q q q q q qq q q q qqq q q q q q q qq q q q q q q q qq q q qq q q q qq qq q q q qq q q q q q q q qqqq q q qq q q q q q qq q qq q q q q q q qq q q q q q q q q q qq q q q qq q q qq q q q qq q q q q q q q q q q qq qq q qq q q q q qq q q q qq q q q q qqq q q q q q q q q q q q qq qq q q 600 700 800 500 q q q qq q 142 140 138 136 134 132 400 q q qq qq 300 q 130 128 66 68 70 72 74 76 200 60 62 64 x (b). Using the plot3d command. 9 y A complete example on kriging using gstat: We will use again the soil data from the Maas river. Here is some background. The actual data set contains many variables but here we will use the x, y coordinates and the concentration of lead and zinc in ppm at each data point. The motivation for this study is to predict the concentration of heavy metals around the banks of the Maas river in the area west of the town Stein in the Netherlands. These heavy metals were accumulated over the years because of river pollution. Here is the area of study: 10 You can access the data at > a <- read.table("http://www.stat.ucla.edu/~nchristo/statistics_c173_c273/ soil.txt", header=TRUE) # Save the original image function: > image.orig <- image To load the gstat package type > library(gstat) First, we will compute the descriptive statistics of the data set, construct the stem-and-leaf plots, histograms, and boxplots: > > > > > > > stem(a$lead) boxplot(a$lead) hist(a$lead) stem(a$zinc) boxplot(a$zinc) hist(a$zinc) summary(a) Transform the data (logarithm transformation): > > > > > > > > log_lead <- log10(a$lead) log_zinc <- log10(a$zinc) stem(log_lead) boxplot(log_lead) hist(log_lead) stem(log_zinc) boxplot(log_zinc) hist(log_zinc) #Create a gstat object; > g <- gstat(id="log_lead", formula = log(lead)~1, locations = ~x+y, data = a) #Plot the variogram: > plot(variogram(g), main="Semivariogram of the log_lead") #Fit a model variogram to the sample variogram: > v.fit <- fit.variogram(variogram(g), vgm(0.5,"Sph",1000,0.1)) > plot(variogram(g),v.fit) #Note: The values above were the initial values for the partial sill, #range, and nugget. Then the function fit.variogram uses a minimization #procedure to fit a model variogram to the sample variogram. Type v.fit 11 #to get the estimates of the model parameters. > v.fit model psill range 1 Nug 0.05156252 0.0000 2 Sph 0.51530678 965.1506 #There are different weights you can use in the minimization procedure. The #default (the one used above) is $N_h/h^2$ where $N_h$ is the number of pairs #and $h$ the separation distance. You can chose the type of weights by using #the argument fit.method=integer, where integer is a number from the table #below: fit.method 1 2 6 7 weights N_h N_h/gamma(h;theta)^2 (Cressie's weights) OLS (no weights) N_h/h^2 (default) #Use kriging to estimate the value of log(lead) at the grid values. #First we create the grid. > x.range <- as.integer(range(a[,1])) > x.range > y.range <- as.integer(range(a[,2])) > y.range > grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=50), y=seq(from=y.range[1], to=y.range[2], by=50)) #We want now to use kriging to predict log(lead) at each point on the grid: > pr_ok <- krige(id="log_lead",log(lead)~1, locations=~x+y, model=v.fit, data=a, newdata=grd) #To find what the object pr_ok contains type: > names(pr_ok) [1] "x" "y" "log_lead.pred" "log_lead.var" #To see the predicted values you type: > pr_ok$log_lead.pred #And the kriging variances: > pr_ok$log_lead.var 12 The plot of the sample variogram: > plot(variogram(g), main="Semivariogram of the log_lead") Semivariogram of the log_lead q 0.6 q q q q 0.5 q q q q q 0.4 semivariance q q 0.3 q 0.2 q 0.1 q 500 1000 1500 distance The fitted spherical variogram to the sample variogram: > plot(variogram(g),v.fit) Fitted spherical semivariogram q 0.6 q q q q 0.5 q q q q q 0.4 semivariance q q 0.3 q 0.2 q 0.1 q 500 1000 1500 distance 13 Here is the grid for the kriging predictions: > plot(grd) Grid for the kriging predictions qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq y 330000 331000 332000 333000 178500 179000 179500 180000 x 180500 181000 Load the library scatterplot3d: > library(scatterplot3d) > scatterplot3d(pr_ok$x, pr_ok$y, pr_ok$log_lead.pred, main="Predicted values") Predicted values q qq qq qq q q q qq q qqq q qqq qq q qqq q qqqq q q q q q qq q q q q q q q qq q q qq q q q qqq qqq q qq q qqq q qq qqq qq q q q qqq qqqq qq qq q qq qqqq qqqq q q q q qqq q qqqq q q q q q q qq q q q q q q qq qq q qq q q q q q qq q q q qqq qq qqqq qq qqq q q q q qqqq q q qq qq q q q q q qq q q q q q q q q q qqq qq qq q q qq q q q qq q qq q q q qq q q q q q qqq q q q q q q qq qq q qqqq q qqqqqqq q q q q qq q q qq q q q q q q q q qq qqq q q qq qq q qq q q qqq q q qq qqq q q q q q q q q qq q q q qq qq q q qq q q qq q qq q q qqqq qq q qq qqqqq qq q q q q qqq q q q q q q q q qq qq qq q q qq q qqqq q q qqq q q qq q qq q q qqq qq qqq qq q q qq q qq qqq qqqqqqqqqq q q qqq qq qqqqq q q qqq q qqq q q q q qqqqqq q q q q q qqqqqq q q q q qqqq q qq q qqq q q q q q q q q q qq q q q q qq q q q qqqq q q q qq q qqqq q q qq q q q q qqq q q q q q qqqqqq q q qq q qqqqq q q q q q qqq q qqq q q qq qq q q q q qq q q q qq q q q qqq q qq q q qq qq q q qqq q q q qq q qqqqqqq qq q q q qqqqqq qq qq q q q q q q qqq q q qq qq qqqqqq qq q qq q qqqqq q q q q q q qq qq qq qq qq q q qqqqqq qq q qqqq qq qqq q qq q qq q qqqqq q q qqq q qq qq q q q q qqqq q q qqq qqq q q q qq qq qq qqqqqq q qq q q q q qq q qq q q qq q qq qqq qq q q q q qqq q q q qq qqqq q q q qq qqq qq q qqqqq q qqq q q q qqq q q q qqqqqqqqq qq qqqq qqq q qq qqqqq q q qq q q qqqqq qqqqqqqqqqqqqqqqqqqqqqqqq qqq qq q qqq q q q q q q qq q q qq qq q qq q qq q q qq q qqqq q qq qqqqqqqqqqqqqqqq qqqqqqqqqqq q q qqq qq qqq qqqqqqqqqqqqqqqqqqqqqqqqqq qqqq q q qq qqqqqqqq q q q q q q qq qqqqq q q qq q q q q qqq qqqq qqqq qqqqqqq q qqq q q qqqqqqqqqqqqqqqqq qqqqq q qq qq qqqqqqqqqqqqqqqqqqqqqqqqqq qqqqqqqqqqqqqqqqq qqqqqqq qq qq q qq q q qqqqqqqqqqqqqq q qqqqqqq q qq qqq qq qq q qqq q q qq q qq q q q q q qq q qq qqqqq qq qqqqqqqqqqqqqqqqqqqqqqq qqq qqqqqqqqq qq qq q q qq q qqqqqqqqqq qqqqq q q qqqqqqqqqqqqqqq q q q q q q q qq qqqqqqqq q qqqqqqqqqqqqqqqqqq qq qqqqqqqqqqqqqqqqqqqqqqqq qq q q qqqq q qqqqqqqqqqqqqqqqqqqq qq q q qq q qq q q qqq qq q qq q qq qqqqqqq q qq q q qq q q qqqqqqqq q qq q q q q q q q q q q qq q qqqqqq q q qq q q q q q qq qqq qqq qqqq qqqqqqqq qqqqq qqq q q q q q q q q qqq q q qqqq qq q q qqqqqqq q q q q qq q q q qqqqqqqqqqq q qqqq q qqqqqqqqq qqqqqqq q q q q q qq qq qq q q q q q qqqqq q qq qqqqqqqqqq q q qq q q q q q q qqqq q q q q qq q q qq q qqq qqqqqq q qq qqq q q qq qqqqqqqq q qq q q q q q q qqqqq q q q qqqqq q qqqq q q q q q q qq q q q q qq qq qq q q q q q q q q q q q q qqq qqq q q q qqq qq q qq q q q q q q qqqq q q q qqqq q qq q q q qqq qq qq q qq q qq q q qqqq q q qqqq q qqq qq q q q q qq q q qqq q q q q qq qq q q q q q q qqq q q q qqq q q qq q q q qq q qq q q q q q q q q qq q qq q qq qqq q q q q qqqq q q qqq q qqqq qq q q q qqqqq q q q q q qq q qqq q qq qqq q q q q qqqqq qq q q q qqq q qqq q qqq q qq q qqq qqqqqq qq q q qqqqqqq q qq q qq qq q q qq q qqqqqqqq qq q q q qq qq q q q q qq qqq qq q q q q q qq qq qqqq q q q q qq q qq q qq q q q q q qq q q q q q q qq qq qq q q q q qq q q q q q q q q q qqqq qqqq qq q qqqqqq q q q q q qq qqqqq q q q q q q q qq qq q q q q qqq q q q q q q q q q qq qq qq q qq qqqq q qqqqq qqqqqqqq q qqqq qqq q qq qq qqqq qq q q q qqqqq q q q q q qq q q q qqqq qqqqqqqqqq q qq q q qq q q q qqqqq q qq qq q qqq q qq q q q qq q q q q qqq q qqq q q qqq qqqq q q q qq q q qqqq qq q qqq q q q q q qq q q q q qqqqq q q qq qqqq q qqqqqqqq qqqqqqq qq qqq q q qq q qqq qq qqq q q qq qqqq q q q qqqqq q qq q qqq q q q q qq qq q qqq q qq qq q q q q q q qq q qqqqq q qq qq qqq q qqq q qq qqq q qqqq q qq q q q q q q q q qqq q qq qq q qq qqqqq q qqqqqqqqqq qqqqqqqqqq qqqqqqqqq q qq qq q qqqqqqq q qq qqqqqq qqq q q q q q q qq q qqqqqqqq q qqqqqqqq qqq q qq q q qqqqqqq qqqqqqqq q qq qqqqqq q q q q q qqq q qq q q q qq q q qqqq qqqqqqq q q q q q q q q q qqq qqq q q qqqqqqqqqqqqqqqqq qqqqqqq q qqq q q qqqq qqqqqq q qqqq q q q q q q q q q qq q qq q q qqq qq q q q q q q q q q q q qq qqq q qq q q qq q qq qqqqqqqqqqqqq q qqqqqq qq q qqqqqq q q qqqqq q q q qqqqq q qqq qq qq q q q q q q q q q qqq q q qqqqqq qqq qqq qqqqqqqqqqq q qqqqq q q qq qq qq qq q qq qqqqq q qqqqqq q q q qqqq q qq q q q q qq q q q qqqq q qq qqqq qq q qqqqqqqqqqq q q qq q q q q qq q qq q q qqqq q qq qq q qqqqqqqqqqqqq qqqq q qqqqqq q q qqqqqqqqqq q q q q q qq q q qq qq q q q q qqqqqqqq qqq q qq q q qq q qq q qqq qqqqq q qqqqqqqqqqqqqqq qqq q qq q qqq q q qq qqq q q q qq q qq q q qq q q q q q q q q qq qqq qq qq q q qq q q qqq qq q qqq qq q qqqqqqqqqq q q q qq qq q q q q q q q q qq q q q qqq q q q qq q q q q qq q q qqq q q q qq q q q q q q qq qqq q q q qqq qq q q q q q qqqq qqqqqqq qqqqqq q q q q qqqqq q qq q qq qq qqq qqq q q qq q q qq q q qq q q q qq qq qq q qq q q qq q qq qqqq q q q q q q qqqqq qq qqq qq q qq q qqqqq qq qq q q q qq q q q q qq q q q qqqq q q q q q q qqq q qq q qq q q q q q q q qq q q q qq q q q q q q qq qq qqq q qq qq q q q q q q q q qq q qqqq qqq qq q q q q q qqq q qqq q qqqqq qqq q q q q qqqqqqq q qqq q q q qq q q qq q q q q q q qq q q qq q q q q q q q q q q qq q qq qq q qqqq q q q q qq q qqqqq qq qqq qq q 334000 qq q q qq q qq q q qq q qq qqqq q qqq q qq q qq qqqqqqqqqq q qq q q qqqqq qqqqqq q qqqqqq q qqqqqq qqqq qq qqqqqqq q qq qq qqqqqqqq q qqqq q q qq q qqqq q qqqqq q q q qqq q qqqqq q qq qq q qqqqqqqqqq q qq qq qq q q q q q q qqq q qq q q q qq q q q q qq q qq qq q qq qq q qq q q qqqq q q q q qq q q q q qq q q qq qq q qqqqqqq q qqqq q qq q q q q qq q q q qq qqq qqqqqq q q q q q q qq qq q q qq 333000 q qqq q qq qqqq qq qq q qqq q q q qqq q q q qq q qq q qqq q q q qq q q qq qqqqqqq qq qqqqq qqqqqqq qqqq qq q q q q qqq q q q q qq q qqq q qq qq q q qqq q q q qqq q qq q qq qq q qq q q qq q q qq q q qq qqq qq qq q q qq 332000 q qq q q pr_ok$log_lead.pred 5.0 5.5 6.0 6.5 4.0 331000 330000 3.5 329000 178500 179000 179500 180000 180500 181000 181500 pr_ok$x 14 pr_ok$y 4.5 Load the library lattice to use the levelplot function: > library(lattice) > levelplot(pr_ok$log_lead.pred~x+y,pr_ok, aspect ="iso", main="Ordinary kriging predictions") Ordinary kriging predictions 333000 6.0 5.5 332000 5.0 y 331000 4.5 4.0 330000 3.5 179000 179500 180000 180500 181000 x > levelplot(pr_ok$log_lead.var~x+y,pr_ok, aspect ="iso", main="Ordinary kriging variances") Ordinary kriging variances 0.6 333000 0.5 0.4 332000 y 0.3 331000 0.2 0.1 330000 0.0 179000 179500 180000 180500 181000 x 15 > levelplot(sqrt(pr_ok$log_lead.var)~x+y,pr_ok, aspect ="iso", main="Ordinary kriging standard errors") Ordinary kriging standard errors 0.8 333000 0.7 0.6 0.5 332000 0.4 y 0.3 331000 0.2 0.1 330000 0.0 179000 179500 180000 180500 181000 x Use the image function: We will use now the function image.orig that was saved at the beginning (before we loaded the gstat library) to construct similar graphs. # Function to convert a vector to a matrix: vec.2.Mx <- function( xvec, nrow, ncol ) { Mx.out <- matrix(0, nrow, ncol ) for(j in 1:ncol) { for(i in 1:nrow) { Mx.out[i, j] <- xvec[ (j-1)*nrow + i ] } } return( Mx.out ) } 16 Now we can use the function to create a matrix (qqq below): > qqq <- vec.2.Mx( pr_ok$log_lead.pred, length(seq(from=x.range[1], to=x.range[2], by=50)), length(seq(from=y.range[1], to=y.range[2], by=50)) ) Much easier the vector of the predicted values can be collapsed into a matrix with the matrix function: #Collapse the predicted values into a matrix: qqq <- matrix(pr_ok$log_lead.pred, length(seq(from=x.range[1], to=x.range[2], by=50)), length(seq(from=y.range[1], to=y.range[2], by=50))) And we can use the original image function to create the raster map of the predicted values: > image.orig(seq(from=x.range[1], to=x.range[2], by=50), seq(from=y.range[1], to=y.range[2], by=50), qqq, xlab="West to East", ylab="South to North", main="Predicted values") > points(a) #The data points can be plotted on the raster map. Predicted values q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q 333000 q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q 332000 q q q q q q q q q South to North q q q q q q 331000 q q q qq q q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q 330000 179000 179500 180000 West to East 180500 181000 17 We can also create the raster map of the kriging variances: > qqq <- vec.2.Mx( pr_ok$log_lead.var, length(seq(from=x.range[1], to=x.range[2], by=50)), length(seq(from=y.range[1], to=y.range[2], by=50)) ) > image.orig(seq(from=x.range[1], to=x.range[2], by=50), seq(from=y.range[1], to=y.range[2], by=50), qqq, xlab="West to East", ylab="South to North", main="Kriging variances") Kriging variances South to North 330000 331000 332000 333000 179000 179500 180000 West to East 180500 181000 18 And finally the raster map of the kriging standard errors: > qqq <- vec.2.Mx( sqrt(pr_ok$log_lead.var), length(seq(from=x.range[1], to=x.range[2], by=50)), length(seq(from=y.range[1], to=y.range[2], by=50)) ) > image.orig(seq(from=x.range[1], to=x.range[2], by=50), seq(from=y.range[1], to=y.range[2], by=50), qqq, xlab="West to East", ylab="South to North", main="Kriging standard errors") Kriging standard errors South to North 330000 331000 332000 333000 179000 179500 180000 West to East 180500 181000 Exercise: Repeat for log(zinc). 19 ...
View Full Document

This note was uploaded on 02/11/2012 for the course STATS c173/c273 taught by Professor Nicolaschristou during the Spring '11 term at UCLA.

Ask a homework question - tutors are online