day2

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library(sp)
library(gstat)

We plot some graphics

#coordinates(meuse)=~x+y
spplot(meuse,"zinc",do.log=T,colorkey=T)

spplot(meuse,"zinc",do.log=T,colorkey=T)

bubble(meuse,"zinc",do.log=T,key.space = "right")

hist(meuse$zinc,xlab="Zinc concentrations", ylab="frequencies",main="Histogram of zinc")

zinc.log<-log(meuse$zinc)
qqnorm(zinc.log,xlab="Zinc (log)", ylab="cumulative frequencies",main="QQplot of zinc log")

#coordinates(meuse)=~x+y
hscat(log(zinc)~1, meuse, (0:9)*100,pch=6)

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

vario

    np       dist     gamma dir.hor dir.ver   id
1   57   79.29244 0.1234479       0       0 var1
2  299  163.97367 0.2162185       0       0 var1
3  419  267.36483 0.3027859       0       0 var1
4  457  372.73542 0.4121448       0       0 var1
5  547  478.47670 0.4634128       0       0 var1
6  533  585.34058 0.5646933       0       0 var1
7  574  693.14526 0.5689683       0       0 var1
8  564  796.18365 0.6186769       0       0 var1
9  589  903.14650 0.6471479       0       0 var1
10 543 1011.29177 0.6915705       0       0 var1
11 500 1117.86235 0.7033984       0       0 var1
12 477 1221.32810 0.6038770       0       0 var1
13 452 1329.16407 0.6517158       0       0 var1
14 457 1437.25620 0.5665318       0       0 var1
15 415 1543.20248 0.5748227       0       0 var1

vario2<-variogram(log(zinc)~1,meuse, cutoff=1600,width=300)
model<-vgm(0.65,"Exp",385,0.0)
plot(vario2, model,main="Exp")

model<-vgm(0.65,"Sph",385,0.0)
plot(vario2, model,main="Sph")

model<-vgm(0.65,"Gau",385,0.0)
plot(vario2, model,main="Gau")

model<-vgm(0.65,"Mat",385,0.0)
plot(vario2, model,main="Mat")

?vgm
vario.fit<-fit.variogram(vario,model=vgm(1,"Sph",900,1))
                         vario.fit

  model      psill    range
1   Nug 0.05066243   0.0000
2   Sph 0.59060780 897.0209

                         plot(vario, vario.fit)

?vgm
vario.fit<-fit.variogram(vario,model=vgm(1,"Sph",900,1))
                         vario.fit

  model      psill    range
1   Nug 0.05066243   0.0000
2   Sph 0.59060780 897.0209

                         plot(vario, vario.fit)

variodir<-variogram(log(zinc)~1, meuse, alpha=c(0,45,90,135))
variodir.fit<-vgm(0.59,"Sph",1200,0.05,anis=c(45,0.4))
plot(variodir,variodir.fit, as.table=TRUE)

We’re going to calculate the kriging with the model previously fitted (vario.fit)

What happens if you don’t use model in kriging

kriging<-krige(log(zinc)~1,meuse,meuse.grid)
kriging$var1.pred[1:10]
spplot(kriging["var1.pred"])

Less go for the directional model

kriging<-krige(log(zinc)~1,meuse,meuse.grid,model=variodir.fit)

[using ordinary kriging]

kriging$var1.pred[1:10]

 [1] 6.656641 6.742468 6.669015 6.565904 6.829509 6.755367 6.642523 6.512539 6.910671
[10] 6.842304

spplot(kriging["var1.pred"])