0

I am in the midway of understanding how gstat package in R tool implements kriging method. I have understood the calculation of empirical semivariogram and fitting semivariogram models. But I did not understand how it implements the matrix inversion to calculate the weights of the kriging estimators. I have a large data set containing 50000 lat-long-precipitaion triplets. Theoretically inversion of a matrix of size 50000x50000 must be done in order to get the weights. While this large matrix takes several GBs of man memory, which is particularly impractical.

My question is that how krige function does all this within a second?

Regards,

Chandan

Chandan
  • 764
  • 2
  • 8
  • 21

1 Answers1

1

You didn't tell what your computing environment is, but I believe it is safe to say that it didn't solve a 50000 points kriging problem in a second. In order to understand what it did, please provide more information, e.g. the commands you used, and the output gstat gave.

Edzer Pebesma
  • 3,814
  • 16
  • 26
  • Thanks for your reply Edzer. I am using a single PC with 4GB of RAM. I am using the following commands prep_data <- read.table("C:/Users/skg/Downloads/data.txt",header=T) attach(prep_data) coordinates(prep_data) = ~x+y prep <- variogram(v ~ x+y,data=prep_data,width=10,cutoff=100) model1.out <- fit.variogram(prep,vgm(70000,"Sph",40,20000)) predpts <- matrix(c(60,190,225,50,110,185),ncol=2,byrow=T) predpts.g <- data.frame(x=predpts[,1],y=predpts[,2]) coordinates(predpts.g) <- ~x+y g <- gstat(NULL,"new.v",v~1,data=prep_data,model=model1.out) three.pred <- predict(g,predpts.g) Regards, Chandan – Chandan Mar 06 '15 at 14:55
  • what output did gstat give after the `predict` step? – Edzer Pebesma Mar 06 '15 at 22:47
  • The output is the value at three prediction point three.pred coordinates new.v.pred new.v.var 1 (60, 190) 1305.29919 34483.78 2 (225, 50) 311.84032 54113.63 3 (110, 185) 65.89723 57739.50 – Chandan Mar 07 '15 at 04:02
  • Basically I want to know what is inversion method implemented in the krige or predict function. – Chandan Mar 07 '15 at 04:39
  • 1
    LDL (Choleski) decomposition; see here: http://stackoverflow.com/questions/19226816/how-can-i-view-the-source-code-for-a-function – Edzer Pebesma Mar 07 '15 at 11:39
  • Can you please tell me from where you found that its using Choleski decomposition. The predict methods gives 15 methods [1] predict.ar* predict.Arima* [3] predict.arima0* predict.glm [5] predict.HoltWinters* predict.lm [7] predict.loess* predict.mlm* [9] predict.nls* predict.poly* [11] predict.ppr* predict.prcomp* [13] predict.princomp* predict.smooth.spline* [15] predict.smooth.spline.fit* predict.StructTS* – Chandan Mar 07 '15 at 16:20
  • Which of these giving the kriging prediction? – Chandan Mar 07 '15 at 16:20
  • 1
    after you load `gstat`, there is a `predict.gstat` too. For kriging, it finally calls https://r-forge.r-project.org/scm/viewvc.php/pkg/src/gls.c?view=markup&root=gstat which then calls LDLsolve() in the matrix routines. – Edzer Pebesma Mar 08 '15 at 14:06