## ----echo=FALSE, message=FALSE, warning=FALSE--------------------------------- library(RGeostats) rm(list=ls()) constant.define("asp",1) ## ----eval=FALSE--------------------------------------------------------------- ## Exdemo_2D_pollution.table = ## read.table("Pollution.dat",header=T,na.strings="-1") ## ----------------------------------------------------------------------------- data(Exdemo_2D_pollution.table) ## ----------------------------------------------------------------------------- data.db <- db.create(Exdemo_2D_pollution.table, flag.grid=FALSE,ndim=2,autoname=F) ## ----------------------------------------------------------------------------- print(data.db) ## ----------------------------------------------------------------------------- print(data.db,flag.stats=TRUE,names=4) ## ----------------------------------------------------------------------------- data.db = db.locate(data.db,"Zn","z",1) ## ----------------------------------------------------------------------------- plot(data.db,pch=21,bg="red",col="black", title="Zn Samples location") ## ----------------------------------------------------------------------------- data.db = db.sel(data.db,Zn<20) ## ----------------------------------------------------------------------------- print(data.db,flag.stats=T,names="Zn") ## ----------------------------------------------------------------------------- data.db = db.add(data.db,Zn.transformed=Zn+Pb+2) ## ----------------------------------------------------------------------------- print(data.db,flag.print=TRUE) ## ----------------------------------------------------------------------------- data.vario = vario.calc(data.db,lag=1,nlag=10) ## ----------------------------------------------------------------------------- plot(data.vario,npairdw=TRUE,npairpt=1,title="Experimental Variogram") ## ----------------------------------------------------------------------------- data.4dir.vario = vario.calc(data.db,lag=1,nlag=10,dir=c(0,45,90,135)) ## ----------------------------------------------------------------------------- plot(data.4dir.vario,title="Directional variograms") ## ----------------------------------------------------------------------------- plot(data.4dir.vario,idir0=2, title="Directional variogram #2") ## ----------------------------------------------------------------------------- data.model = model.auto(data.vario, struct=c("Spherical","Exponential"), title="Modelling omni-directional variogram") ## ----------------------------------------------------------------------------- data.model ## ----------------------------------------------------------------------------- Exdemo_2D_unique.neigh = neigh.create(ndim=2, type=0) ## ----------------------------------------------------------------------------- Exdemo_2D_unique.neigh ## ----------------------------------------------------------------------------- Exdemo_2D_moving.neigh = neigh.create(ndim=2,nmaxi=10,radius=20, flag.sector=TRUE,nsect=8,nsmax=2) ## ----------------------------------------------------------------------------- Exdemo_2D_moving.neigh ## ----------------------------------------------------------------------------- data.db = xvalid(data.db,data.model, Exdemo_2D_moving.neigh) ## ----------------------------------------------------------------------------- hist(db.extract(data.db,"Xvalid.Zn.*esterr"), nclass=30, main="Histogram of Zn Error", xlab="Cross-validation error", col="blue") ## ----------------------------------------------------------------------------- grid.db = db.grid.init(data.db,nodes=c(100,90)) ## ----------------------------------------------------------------------------- grid.db ## ----------------------------------------------------------------------------- data.db <- db.locerase(data.db,"z") data.db <- db.locate(data.db,"Zn","z") data.db ## ----------------------------------------------------------------------------- grid.db = kriging(data.db,grid.db,data.model, Exdemo_2D_unique.neigh,radix="KU.Part") ## ----------------------------------------------------------------------------- plot(grid.db,name.image="KU.Part.Zn.estim", title = "Estimation - Data subset (Unique Neighborhood)") plot(grid.db,name.contour="KU.Part.Zn.stdev",nlevels=10,add=TRUE) plot(data.db,pch=21,bg="yellow",col="black",add=TRUE) ## ----------------------------------------------------------------------------- plot(grid.db,name.persp="KU.Part.Zn.estim",theta=45,phi=30, zlim=c(3,25), title = "Estimation - Data subset (Unique Neighborhood)") ## ----------------------------------------------------------------------------- data.db = db.sel(data.db) ## ----------------------------------------------------------------------------- grid.db = kriging(data.db,grid.db,data.model, Exdemo_2D_unique.neigh,radix="KU.All") ## ----------------------------------------------------------------------------- plot(grid.db,name.image="KU.All.Zn.estim", title = "Estimation - All Data (Unique Neighborhood)") plot(grid.db,name.contour="KU.All.Zn.stdev",nlevels=10,add=TRUE) plot(data.db,pch=21,cex1=0.5,bg="yellow",col="black",add=TRUE) ## ----------------------------------------------------------------------------- plot(grid.db,name.persp="KU.All.Zn.estim",theta=45,phi=30, zlim=c(3,25), title = "Estimation - All Data (Unique Neighborhood)") ## ----------------------------------------------------------------------------- grid.db = kriging(data.db,grid.db,data.model, Exdemo_2D_moving.neigh,radix="KM.All") ## ----------------------------------------------------------------------------- plot(grid.db,name.image="KM.All.Zn.estim", title = "Estimation - Data subset (Moving Neighborhood)") plot(grid.db,name.contour="KM.All.Zn.stdev",nlevels=10,add=TRUE) plot(data.db,pch=21,cex1=0.5,bg="yellow",col="black",add=TRUE) ## ----------------------------------------------------------------------------- grid.db = neigh.test(data.db,grid.db,data.model, Exdemo_2D_moving.neigh,radix="Moving") ## ----------------------------------------------------------------------------- plot(grid.db,name="Moving.Zn.max.radius", title="Largest Distance") plot(grid.db,name="Moving.Zn.min.dist", title="Smallest Distance") plot(grid.db,name="Moving.Zn.nb.samples", title="Number of selected samples") plot(grid.db,name="Moving.Zn.nb.nonempty.sect", title="Number of non-empty consecutive sectors") ## ----------------------------------------------------------------------------- data.anam = anam.fit(data.db,"Zn",title="Anamorphosis Function") ## ----------------------------------------------------------------------------- print(data.anam) ## ----------------------------------------------------------------------------- data.db = anam.z2y(data.db,"Zn",anam=data.anam) ## ----------------------------------------------------------------------------- print(data.db,flag.stats=TRUE,names="Gaussian.Zn") ## ----------------------------------------------------------------------------- data.g.vario = vario.calc(data.db,nlag=10,lag=1) ## ----------------------------------------------------------------------------- plot(data.g.vario,npairdw=TRUE,npairpt=1, title="Variogram of Zn Gaussian") ## ----------------------------------------------------------------------------- data.g.model = model.auto(data.g.vario,struct=c("Exponential"), title="Model for Zn Gaussian") ## ----------------------------------------------------------------------------- grid.db = simtub(data.db,grid.db,data.g.model, Exdemo_2D_unique.neigh,nbsimu=10,nbtuba=100) ## ----------------------------------------------------------------------------- grid.db = anam.y2z(grid.db,names="Simu.Gaussian.Zn*",anam=data.anam) ## ----------------------------------------------------------------------------- plot(grid.db,name.image="Raw.Simu.Gaussian.Zn.S1", zlim=c(3,13),title="Simulation #1") plot(data.db,pch=21,cex1=0.5,bg.in=1,add=TRUE) plot(grid.db,name.image="Raw.Simu.Gaussian.Zn.S2", zlim=c(3,13),title="Simulation #2") plot(data.db,pch=21,cex1=0.5,bg.in=1,add=TRUE) ## ----------------------------------------------------------------------------- grid.db = db.compare(grid.db,names="Raw.Simu.Gaussian.Zn*",fun="mean") ## ----------------------------------------------------------------------------- plot(grid.db,zlim=c(3,13),title="Mean of Simulations") plot(data.db,pch=21,cex1=0.5,bg.in=1,add=TRUE) ## ----------------------------------------------------------------------------- dummy = correlation(grid.db,"KU.All.Zn.estim","mean", flag.iso=TRUE,flag.diag=TRUE, xlab="Kriging Estimate", ylab="Mean of Simulations", title="Comparing Simulations and Estimation") ## ----------------------------------------------------------------------------- grid.db = db.compare(grid.db,names="Raw.Simu.Gaussian.Zn*",fun="stdv") ## ----------------------------------------------------------------------------- plot(grid.db,name="stdv",title="Dispersion of Simulations") plot(data.db,pch=21,cex1=0.5,bg.in=1,add=TRUE)