Dear Francky,
here is a short script using a data set of RGeostats to explain how to perform universal kriging.
Load the data:
- Code: Select all
data(Exdemo_ext.data)
dat=Exdemo_ext.data #to do shorter
Perform the variography and the modeling:
- Code: Select all
vario=vario.calc(dat)
model1=model.auto(vario,struct=c("K-Bessel","Nugget Effect"))
Create the output grid and the neighborhood (unique):
- Code: Select all
grid=db.create(nx=c(200,200))
unique=neigh.init(type=0)
Perform ordinary kriging:
- Code: Select all
grid=kriging(dat,grid,neigh=unique,model=model1)
plot(grid,pos.legend=1,zlim=c(-4,5))
plot(dat,add=T)
Perform universal kriging with a second order polynomial drift:
- Code: Select all
x11()
grid=kriging(dat,grid,neigh=unique,model=model1,uc=c("1","x","y","xy","x2","y2"))
plot(grid,pos.legend=1,zlim=c(-4,5))
plot(dat,add=T)
It seems that there is not a great difference (it is probably due to the
fact that drift is useless in this specific case and the data control the result.
However we can see some differences far from the data points:
- Code: Select all
grid=db.add(grid,grid[,6]-grid[,4],loctype="z")
plot(grid,pos.legend=1)
plot(dat,add=T)
I hope it helped you.
I am not really aware about universal kriging
Regards,
Nicolas