Dear Jui-Han, I believe the following script does what you are looking after. The idea is to start with the covariance matrix purely defined by the variogram model (i.e. with the joint sill on the diagonal). Rescale this covariance matrix first to correlations and then back to covariances using the marginal variances obtained from the kriging predictor. I will correct the krige0 function in gstat accordingly.
Let me know if this is what you have been looking for. Best wishes, Ben #--R-Script--# # get some data library(gstat) library(sp) data(meuse) coordinates(meuse) = ~x+y data(meuse.grid) gridded(meuse.grid) = ~x+y m <- vgm(0.59,"Sph", 897, 0.05) # do the prediction to get the marginal kriging # variance for each grid cell pred <- krige(log(zinc)~1, meuse, meuse.grid, m) # calculate the covariacne strucutre of the desired # grid: 3103 x 3103 matrix covMat <- variogramLine(m, dist_vector = spDists(meuse.grid, meuse.grid), covariance = TRUE) # rescale to correlation matrix corMat <- cov2cor(covMat) # scale back to the desired variances on the # diagonal, just the same way backwards as cov2cor # does it: covMat <- corMat*matrix(sqrt(pred$var1.var) %x% sqrt(pred$var1.var),3103,3103) covMat[1:5,1:5] # check for correlation cov2cor(covMat[1:5,1:5]) # check prediction variance: diag(covMat) - pred$var1.var #--R-Script--# On 12.12.2013 22:54, Edzer Pebesma wrote: > > > On 12/12/2013 08:20 PM, Jui-Han Chang wrote: >> Hi all, >> >> I am trying to get covariance matrix for kriging predictions. I used krige0 >> function in the gstat library with options fullCovariance=TRUE and >> computeVar=TRUE. I got a covariance matrix but the values seem to be too >> large. I converted the covariance matrix to correlation matrix and got >> values larger than 1. Following are the codes to derive the matrices: >> >> library(gstat) >> data(meuse) >> coordinates(meuse) = ~x+y >> data(meuse.grid) >> gridded(meuse.grid) = ~x+y >> >> m <- vgm(.59, "Sph", 874, .04) >> x <- krige0(log(zinc)~1, meuse, meuse.grid, model = >> m,fullCovariance=TRUE,computeVar=TRUE) >> >> x$var[1:5,1:5] >> >> [,1] [,2] [,3] [,4] [,5] >> 1 0.3108421 0.2794621 0.2883862 0.3019975 0.2528104 >> 2 0.2794621 0.2414045 0.2546174 0.2727450 0.2074984 >> 3 0.2883862 0.2546174 0.2631935 0.2769452 0.2263174 >> 4 0.3019975 0.2727450 0.2769452 0.2863739 0.2499542 >> 5 0.2528104 0.2074984 0.2263174 0.2499542 0.1648699 >> >> cov2cor(x$var[1:5,1:5]) >> >> [,1] [,2] [,3] [,4] [,5] >> 1 1.000000 1.020188 1.008247 1.012201 1.116746 >> 2 1.020188 1.000000 1.010131 1.037331 1.040091 >> 3 1.008247 1.010131 1.000000 1.008764 1.086450 >> 4 1.012201 1.037331 1.008764 1.000000 1.150331 >> 5 1.116746 1.040091 1.086450 1.150331 1.000000 >> >> I am not sure what I am getting here now. Is this covariance matrix for >> kriging predictions? If not does anyone know how to get a covariance >> matrix? Any help will be appreciated. > > Thanks, this is clearly wrong. My feeling is that c0 in the computation > of the kriging predictions covariance matrix needs be specified > differently, not as a constant, however I haven't sorted out how, so far. > > > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Benedikt Gräler ifgi - Institute for Geoinformatics University of Muenster http://ifgi.uni-muenster.de/graeler Phone: +49 251 83-33082 Mail: ben.grae...@uni-muenster.de
signature.asc
Description: OpenPGP digital signature
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
