# NOT RUN {
library(sp)
library(gstat)
# Meuse data:
demo(meuse, echo=FALSE)
# fit a model:
omm <- fit.gstatModel(observations = meuse, formulaString = om~dist,
family = gaussian(log), covariates = meuse.grid)
str(omm@vgmModel)
# write the regression matrix to GeoEAS:
meuse$log_om <- log1p(meuse$om)
write.data(obj=meuse, covariates=meuse.grid["dist"],
outfile="meuse.eas", methodid="log_om")
writeGDAL(meuse.grid["dist"], "dist.rst", drivername="RST", mvFlag="-99999")
makeGstatCmd(log_om~dist, vgmModel=omm@vgmModel,
outfile="meuse_om_sims.cmd", easfile="meuse.eas",
nsim=50, nmin=20, nmax=40, radius=1500)
# compare the processing times:
system.time(system("gstat meuse_om_sims.cmd"))
vgmModel = omm@vgmModel
class(vgmModel) <- c("variogramModel", "data.frame")
system.time(om.rk <- krige(log_om~dist, meuse[!is.na(meuse$log_om),],
meuse.grid, nmin=20, nmax=40, model=vgmModel, nsim=50))
# }
Run the code above in your browser using DataLab