12  meuse dataset

library(sp)
library(gstat)
library(sf)
library(mapview)

data(meuse)
data(meuse.grid)
head(meuse)

       x      y cadmium copper lead zinc  elev       dist   om ffreq soil lime
1 181072 333611    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1
2 181025 333558     8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1
3 181165 333537     6.5     68  199  640 7.800 0.10302900 13.0     1    1    1
4 181298 333484     2.6     81  116  257 7.655 0.19009400  8.0     1    2    0
5 181307 333330     2.8     48  117  269 7.480 0.27709000  8.7     1    2    0
6 181390 333260     3.0     61  137  281 7.791 0.36406700  7.8     1    2    0
  landuse dist.m
1      Ah     50
2      Ah     30
3      Ah    150
4      Ga    270
5      Ah    380
6      Ga    470

head(meuse.grid)

       x      y part.a part.b      dist soil ffreq
1 181180 333740      1      0 0.0000000    1     1
2 181140 333700      1      0 0.0000000    1     1
3 181180 333700      1      0 0.0122243    1     1
4 181220 333700      1      0 0.0434678    1     1
5 181100 333660      1      0 0.0000000    1     1
6 181140 333660      1      0 0.0122243    1     1

meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
meuse

Simple feature collection with 155 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 178605 ymin: 329714 xmax: 181390 ymax: 333611
Projected CRS: Amersfoort / RD New
First 10 features:
   cadmium copper lead zinc  elev       dist   om ffreq soil lime landuse
1     11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1      Ah
2      8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1      Ah
3      6.5     68  199  640 7.800 0.10302900 13.0     1    1    1      Ah
4      2.6     81  116  257 7.655 0.19009400  8.0     1    2    0      Ga
5      2.8     48  117  269 7.480 0.27709000  8.7     1    2    0      Ah
6      3.0     61  137  281 7.791 0.36406700  7.8     1    2    0      Ga
7      3.2     31  132  346 8.217 0.19009400  9.2     1    2    0      Ah
8      2.8     29  150  406 8.490 0.09215160  9.5     1    1    0      Ab
9      2.4     37  133  347 8.668 0.18461400 10.6     1    1    0      Ab
10     1.6     24   80  183 9.049 0.30970200  6.3     1    2    0       W
   dist.m              geometry
1      50 POINT (181072 333611)
2      30 POINT (181025 333558)
3     150 POINT (181165 333537)
4     270 POINT (181298 333484)
5     380 POINT (181307 333330)
6     470 POINT (181390 333260)
7     240 POINT (181165 333370)
8     120 POINT (181027 333363)
9     240 POINT (181060 333231)
10    420 POINT (181232 333168)

mapview(meuse, zcol = "zinc",  map.types = "CartoDB.Voyager")

meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"),
                       crs = 28992)
mapview(meuse.grid,  map.types = "CartoDB.Voyager")

12.1 Variogram Cloud

# Load packages
library(sp)
library(gstat)

# Load Meuse data
data(meuse)
coordinates(meuse) <- ~x+y

# Log-transform zinc
meuse$logZinc <- log(meuse$zinc)
# Compute variogram cloud (all pairs, raw semivariance)
cloud <- variogram(logZinc~1, meuse, cloud = TRUE)

# View first few points
head(cloud)

       dist       gamma dir.hor dir.ver   id left right
1  70.83784 0.006065804       0       0 var1    2     1
2 118.84864 0.109534743       0       0 var1    3     1
3 141.56624 0.167153095       0       0 var1    3     2
4 259.23927 0.952808244       0       0 var1    4     1
5 282.85155 1.110920725       0       0 var1    4     2
6 143.17123 0.416229664       0       0 var1    4     3

# Plot variogram cloud
plot(cloud$dist, cloud$gamma,
     xlab = "Distance (m)",
     ylab = "Semivariance",
     main = "Variogram Cloud",
     pch = 19, col = "blue")

12.2 h-scatterplot

library(sp)
data(meuse)
coordinates(meuse) = ~x+y
hscat(log(zinc)~1, meuse, c(0, 80, 120, 250, 500, 1000))

12.3 Sample variogram

v <- variogram(log(zinc) ~ 1, data = meuse)
plot(v)