library(gstat)
#gstat’ version 1.0-21
library(sp)
# meuse est desormais meuse.all
data(meuse.all)
meuse<-meuse.all
class(meuse)
## [1] "data.frame"
names(meuse)
##  [1] "sample"      "x"           "y"           "cadmium"     "copper"     
##  [6] "lead"        "zinc"        "elev"        "dist.m"      "om"         
## [11] "ffreq"       "soil"        "lime"        "landuse"     "in.pit"     
## [16] "in.meuse155" "in.BMcD"
coordinates(meuse) = ~x+y
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
summary(meuse)
## Object of class SpatialPointsDataFrame
## Coordinates:
##      min    max
## x 178605 181390
## y 329714 333611
## Is projected: NA 
## proj4string : [NA]
## Number of points: 164
## Data attributes:
##      sample          cadmium           copper            lead       
##  Min.   :  1.00   Min.   : 0.000   Min.   : 14.00   Min.   : 27.00  
##  1st Qu.: 41.75   1st Qu.: 0.800   1st Qu.: 23.00   1st Qu.: 68.75  
##  Median : 82.50   Median : 1.900   Median : 29.50   Median :116.00  
##  Mean   : 82.50   Mean   : 3.109   Mean   : 39.42   Mean   :148.55  
##  3rd Qu.:123.25   3rd Qu.: 3.725   3rd Qu.: 48.00   3rd Qu.:201.75  
##  Max.   :164.00   Max.   :18.100   Max.   :128.00   Max.   :654.00  
##                                                                     
##       zinc             elev            dist.m             om        
##  Min.   : 107.0   Min.   : 0.000   Min.   :  10.0   Min.   : 1.000  
##  1st Qu.: 191.8   1st Qu.: 7.390   1st Qu.:  80.0   1st Qu.: 5.000  
##  Median : 307.5   Median : 8.124   Median : 270.0   Median : 6.550  
##  Mean   : 464.6   Mean   : 7.775   Mean   : 294.2   Mean   : 7.291  
##  3rd Qu.: 662.5   3rd Qu.: 8.915   3rd Qu.: 450.0   3rd Qu.: 8.950  
##  Max.   :1839.0   Max.   :10.520   Max.   :1000.0   Max.   :17.000  
##                                                     NA's   :2       
##      ffreq            soil            lime           landuse  
##  Min.   :1.000   Min.   :1.000   Min.   :0.0000   W      :54  
##  1st Qu.:1.000   1st Qu.:1.000   1st Qu.:0.0000   Ah     :42  
##  Median :1.000   Median :1.000   Median :0.0000   Am     :22  
##  Mean   :1.604   Mean   :1.463   Mean   :0.2988   Fw     :10  
##  3rd Qu.:2.000   3rd Qu.:2.000   3rd Qu.:1.0000   Ab     : 8  
##  Max.   :3.000   Max.   :3.000   Max.   :1.0000   (Other):27  
##                                                   NA's   : 1  
##    in.pit        in.meuse155      in.BMcD       
##  Mode :logical   Mode :logical   Mode :logical  
##  FALSE:156       FALSE:9         FALSE:66       
##  TRUE :8         TRUE :155       TRUE :98       
##  NA's :0         NA's :0         NA's :0        
##                                                 
##                                                 
## 
coordinates(meuse)[1:5,]
##        x      y
## 1 181072 333611
## 2 181025 333558
## 3 181165 333537
## 4 181298 333484
## 5 181307 333330
bubble(meuse, "zinc", 
       col=c("#00ff0088", "#00ff0088"), main = "zinc concentrations (ppm)")

data(meuse.grid)
summary(meuse.grid)
##        x                y              part.a           part.b      
##  Min.   :178460   Min.   :329620   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:179420   1st Qu.:330460   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :179980   Median :331220   Median :0.0000   Median :1.0000  
##  Mean   :179985   Mean   :331348   Mean   :0.3986   Mean   :0.6014  
##  3rd Qu.:180580   3rd Qu.:332140   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :181540   Max.   :333740   Max.   :1.0000   Max.   :1.0000  
##       dist        soil     ffreq   
##  Min.   :0.0000   1:1665   1: 779  
##  1st Qu.:0.1193   2:1084   2:1335  
##  Median :0.2715   3: 354   3: 989  
##  Mean   :0.2971                    
##  3rd Qu.:0.4402                    
##  Max.   :0.9926
class(meuse.grid)
## [1] "data.frame"
coordinates(meuse.grid) = ~x+y
class(meuse.grid)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
gridded(meuse.grid) = TRUE

plot(log(zinc)~sqrt(dist.m), meuse)
abline(lm(log(zinc)~sqrt(dist.m), meuse))

#Variogram -> 
lzn.vgm = variogram(log(zinc)~1, meuse)
lzn.vgm
##     np      dist     gamma dir.hor dir.ver   id
## 1   60   80.0948 0.1877905       0       0 var1
## 2  336  164.1146 0.2730703       0       0 var1
## 3  461  267.7307 0.3481733       0       0 var1
## 4  529  372.8053 0.4576170       0       0 var1
## 5  623  478.5473 0.4885060       0       0 var1
## 6  631  585.6598 0.5858425       0       0 var1
## 7  664  692.5486 0.5763304       0       0 var1
## 8  660  795.9618 0.6186735       0       0 var1
## 9  677  902.8950 0.6515299       0       0 var1
## 10 625 1011.8034 0.6692082       0       0 var1
## 11 567 1117.9600 0.7190128       0       0 var1
## 12 542 1221.2978 0.6190030       0       0 var1
## 13 518 1328.9376 0.6803998       0       0 var1
## 14 514 1436.9733 0.6031344       0       0 var1
## 15 462 1542.2698 0.6053041       0       0 var1
lzn.fit = fit.variogram(lzn.vgm, model = vgm(1, "Sph", 900, 1))
lzn.fit
##   model     psill    range
## 1   Nug 0.1229497   0.0000
## 2   Sph 0.5232415 880.8991
plot(lzn.vgm, lzn.fit)

# Kriging
lzn.kriged = krige(log(zinc)~1, meuse, meuse.grid, model = lzn.fit)
## [using ordinary kriging]
spplot(lzn.kriged["var1.pred"])