Unformatted text preview:

University of California, Los AngelesDepartment of StatisticsStatistics C173/C273 Instructor: Nicolas ChristouOrdinary kriging usinggeoRandgstatIn this document we will discuss kriging using the R packages geoR and gstat. We will usethe numerical example from last lecture. Here it is:A simple example:Consider the following datasix y z(si)s161 139 477s263 140 696s364 129 227s468 128 646s571 140 606s673 141 791s775 128 783s065 137 ???Our goal is to estimate the unknown value at location s0. Here is the x − y plot:●●●●●●●●62 64 66 68 70 72 74128 130 132 134 136 138 140x coordinatey coordinates0s1s2s3s4s5s6s71For these data, let’s assume that we use the exponential semivariogram model with param-eters c0= 0, c1= 10, α = 3.33.γ(h) = c0+ c1(1 − e−hα) = 10(1 − e−h3.33).which is equivalent to the covariance functionC(h) =(c0+ c1, h = 0c1e−hα, h > 0⇒ C(h) =(10, h = 010e−h3.33, h > 0The predicted value at location s0is equal to:ˆz(s0) =nXi=1wiz(si) = 0.174(477) + · · · + 0.086(783) = 592.59.And the variance:σ2e=nXi=1wiγ(si− s0) + λ = 0.174(7.384) + · · · + 0.086(9.823) + 0.906 = 8.96.Kriging usinggeoR:We will use now the geoR package to find the same result. First we read our data as ageodata 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 geodatacov.model The model we are usingcov.pars The parameters of the model (partial sill and range)nugget The value of the nugget effectlocations The coordinates (x, y) of the points to be predictedThe 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.9602942Suppose now we want to predict the value at many locations. The following commands willproduce 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 thatpredicts 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 210ksline: kriging location: 101 out of 210ksline: kriging location: 201 out of 210ksline: kriging location: 210 out of 210Kriging performed using global neighbourhoodWe 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.2804173To construct the raster map we type:image(q, val=q$predict)Or simply:image(q)Here is the plot:65 70 75128 130 132 134 136 138 140X CoordY CoordAnd here is the plot with the data points:points(a)65 70 75128 130 132 134 136 138 140X CoordY Coord●●●●●●●4We can construct a raster map of the variances:> image(q, val=q$krige.var)Here is the plot:65 70 75128 130 132 134 136 138 140X CoordY CoordAnd also we can construct a raster map of the standard errors:> image(q, val=sqrt(q$krige.var))Here is the plot:65 70 75128 130 132 134 136 138 140X CoordY Coord5The 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")x coordinateycoordinatePredictedvalues ofzPerspective plot of the predicted valuesAnd 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")x coordinateycoordinatePredicted values of zPerspective plot of the standard errors6Commnets:The argument locations for the function ksline can be a matrix or a data frame. So farwe 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 usedas follows:> kc <- krige.conv(b, loc=in_mat,krige=krige.control(cov.pars=c(10, 3.33), nugget=0))krige.conv: model with constant meankrige.conv: Kriging performed using global neighbourhoodOr 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 meankrige.conv: Kriging performed using global neighbourhoodIn case you have a variofit output you can use it as an input of the argument krige asfollows (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))7Kriging usinggstat:We will use now the gstat package to find the same result. First we read our data andcreate 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 usedas data frames,


View Full Document

UCLA STATS C173 - Lecture Notes

Download Lecture Notes
Our administrator received your request to download this document. We will send you the file to your email shortly.
Loading Unlocking...
Login

Join to view Lecture Notes and access 3M+ class-specific study document.

or
We will never post anything without your permission.
Don't have an account?
Sign Up

Join to view Lecture Notes 2 2 and access 3M+ class-specific study document.

or

By creating an account you agree to our Privacy Policy and Terms Of Use

Already a member?