[R-sig-Geo] Ordinary co-kriging with gstat - computing time/syntax error?
Edzer Pebesma
edzer.pebesma at uni-muenster.de
Mon Apr 10 19:03:35 CEST 2017
You need to set nmax for each variable, as a parameter of the gstat()
function call. I don't see nmax mentioned anywhere in the documentation
of gstat::predict.
On 10/04/17 18:16, Mercedes Román wrote:
> Dear R-sig-Geo members,
>
> I want to perform co-kriging on the residuals on two regression models with
> gstat. I am following some tutorials (Rossiter, 2012; Malone et al., 2017),
> but I may have a mistake in the script. The computing time seems too long
> (I know it takes a lot of calculation time, but in four days it never
> finished), even when I set nmax=20. Originally, I wanted to perform a test
> on ~ 11000 points, but I am trying now with a subset of 200 points. My
> training dataset has ~ 27000 observations. I would like to produce a
> reproducible example, but perhaps just the sintaxis of the script may be
> enough to have a clue of my problem. My computer has an Intel®Xeon® CPU
> E5-2609 v2 processor, and 32Go RAM memory.
>
> Let’s say my two response variables are X1 and X2, and the dataframes with
> the coordinates and covariates are “training” and “test”.
>
> ### Calculate the predictions on the training dataset (predictors is a
> dataframe with the covariates)
> training$pred.X1 <- predict(modelX1, newdata = predictors)
> training$pred.X2 <- predict(modelX2, newdata = predictors)
> ### Calculate residuals
> training$res.X1<- training$X1 - training$ pred.X1
> training$res.X2 <- training$X2 - training$pred.X2
> ### Transform to spatial
> training.sp <- training
> coordinates(training.sp) <- ~ x + y
> proj4string(training.sp) <- CRS("+init=epsg:2154")
>
> ### Define gstat object and compute experimental semivariograms for X1 and
> X2 (for visual examination)
> X1.g <- gstat(formula = res.X1~1, data = training.sp)
> X1.svar <- variogram(X1.g, width = 5000, cutoff=250000)
> X1.ivgm <- vgm(nugget=0.56, range=14229, psill=0.65, kappa=10, model="Mat")
> ## initial variogram model
> X1.vgm <- fit.variogram(X1.svar, model=X1.ivgm, fit.method=7)
> plot(X1.svar, X1.vgm)
>
> X2.g <- gstat(formula = res.X2~1, data = training.sp)
> X2.svar <- variogram(X2.g, width = 5000, cutoff=250000)
> X2.ivgm <- vgm(nugget=0.19, range=13954, psill=0.21,kappa=10, model="Mat")
> X2.vgm <- fit.variogram(X2.svar, model=X2.ivgm, fit.method=7)
> plot(X2.svar, X2.vgm)
>
> ### Define gstat object to compute the direct variograms and covariogram
> g <- gstat( NULL, id="res.X1", formula = res.X1~1, data = training.sp)
> g <- gstat( g, id="res.X2", formula = res.X2~1, data = training.sp)
> v.cross <- variogram(g, width=5000, cutoff = 250000)
> #### Fill parameters
> g <- gstat(g, id = "res.X2", model = X2.vgm, fill.all=T)
> ### fit LMCR
> g <- fit.lmc(v.cross, g)
> ### Predict at the test locations (test.sp only indicates the coordinates)
> test.sp <- test
> coordinates(ensemble.sp) <- ~ x + y
> proj4string(ensemble.sp) <- CRS("+init=epsg:2154")
> k.c_residuals <- predict(g, test.sp, nmax=20)
>
> I am wondering if there is something missing in the predict command… I do
> not need to include the training data in it?
> I have been reading previous posts from R-sig-Geo, but I still did not find
> the solution. I would appreciate any advice you could give me.
>
> Thanks,
>
> Mercedes Roman
>
> INRA
>
> Unité de Service InfoSol
>
> 2163, avenue de la Pomme de Pin
>
> CS 40001 Ardon
>
> 45075 ORLEANS cedex 2
>
> France
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software: http://www.jstatsoft.org/
Computers & Geosciences: http://elsevier.com/locate/cageo/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20170410/28442d14/attachment.sig>
More information about the R-sig-Geo
mailing list