This article shows how to krig soil data with the function krig() developed by Graham Zemunik (see ?krig()).

Compare parameters

Krige with automated parameters.

auto <- krig(soil, var = "m3al", quiet = TRUE)
#> Guessing: plotdim = c(1000, 500)
summary(auto)
#> var: m3al 
#> df
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1250 obs. of  3 variables:
#>  $ x: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ y: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z: num  950 944 939 933 928 ...
#> 
#> df.poly
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1250 obs. of  3 variables:
#>  $ gx: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ gy: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z : num  949 943 938 933 927 ...
#> 
#> lambda
#> 'numeric'
#>  num 1
#> 
#> vg
#> 'variogram'
#> List of 20
#>  $ u               : num [1:12] 30.3 42.9 51.1 60.9 86.5 ...
#>  $ v               : num [1:12] 60522 54239 174166 55226 44193 ...
#>  $ n               : num [1:12] 21 25 5 81 94 57 155 167 247 288 ...
#>  $ sd              : num [1:12] 96683 58618 214627 60326 52526 ...
#>  $ bins.lim        : num [1:31] 1.00e-12 2.00 2.38 2.84 3.38 ...
#>  $ ind.bin         : logi [1:30] FALSE FALSE FALSE FALSE FALSE FALSE ...
#>  $ var.mark        : num 55039
#>  $ beta.ols        : num 2.62e-08
#>  $ output.type     : chr "bin"
#>  $ max.dist        : num 320
#>  $ estimator.type  : chr "classical"
#>  $ n.data          : int 100
#>  $ lambda          : num 1
#>  $ trend           : chr "cte"
#>  $ pairs.min       : num 5
#>  $ nugget.tolerance: num 1e-12
#>  $ direction       : chr "omnidirectional"
#>  $ tolerance       : chr "none"
#>  $ uvec            : num [1:30] 1 2.19 2.61 3.11 3.7 ...
#>  $ call            : language variog(geodata = geodata, breaks = breaks, trend = trend, pairs.min = 5)
#> 
#> vm
#> 'variomodel', variofit'
#> List of 17
#>  $ nugget               : num 56634
#>  $ cov.pars             : num [1:2] 87783 2105
#>  $ cov.model            : chr "gaussian"
#>  $ kappa                : num 0.5
#>  $ value                : num 1.49e+11
#>  $ trend                : chr "cte"
#>  $ beta.ols             : num 2.62e-08
#>  $ practicalRange       : num 3644
#>  $ max.dist             : num 320
#>  $ minimisation.function: chr "optim"
#>  $ weights              : chr "npairs"
#>  $ method               : chr "WLS"
#>  $ fix.nugget           : logi FALSE
#>  $ fix.kappa            : logi TRUE
#>  $ lambda               : num 1
#>  $ message              : chr "optim convergence code: 0"
#>  $ call                 : language variofit(vario = vg, ini.cov.pars = c(initialVal, startRange), cov.model = varModels[i],      nugget = initialVal)

Try also:

View(auto)

With custom parameters (arbitrary but based on automated kriging parameters).

Compare.

title_auto <- "Automated parameters"
ggplot(as_tibble(auto), aes(x = x, y = y, fill = z)) +
  geom_tile() + 
  coord_equal() +
  labs(title = title_auto)

title_custom <- "Custom parameters"
ggplot(as_tibble(custom), aes(x = x, y = y, fill = z)) +
  geom_tile() + 
  coord_equal() +
  labs(title = title_custom)

Iterate over multiple soil variables

Scale soil data to make their values fall within the same range. This is not crucial, but it will ease visual comparison.

soil_vars <- c("mg", "c", "p")
soil_scaled <- modify_at(soil_fake, .at = soil_vars, scale)

krig_list <- krig(soil_scaled, soil_vars, quiet = TRUE)
#> Guessing: plotdim = c(1000, 460)

summary(krig_list)
#> var: mg 
#> df
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1150 obs. of  3 variables:
#>  $ x: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ y: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z: num  -0.3186 -0.2256 -0.1372 -0.0535 0.0256 ...
#> 
#> df.poly
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1150 obs. of  3 variables:
#>  $ gx: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ gy: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z : num  -0.3186 -0.2256 -0.1372 -0.0535 0.0256 ...
#> 
#> lambda
#> 'numeric'
#>  num 1
#> 
#> vg
#> 'variogram'
#> List of 20
#>  $ u               : num [1:9] 60.9 86.5 103 122.7 146.1 ...
#>  $ v               : num [1:9] 0.858 0.852 1.19 0.856 1.229 ...
#>  $ n               : num [1:9] 7 9 10 10 18 19 36 34 38
#>  $ sd              : num [1:9] 1.204 0.981 1.444 0.835 0.996 ...
#>  $ bins.lim        : num [1:31] 1.00e-12 2.00 2.38 2.84 3.38 ...
#>  $ ind.bin         : logi [1:30] FALSE FALSE FALSE FALSE FALSE FALSE ...
#>  $ var.mark        : num 0.819
#>  $ beta.ols        : num -5.08e-09
#>  $ output.type     : chr "bin"
#>  $ max.dist        : num 320
#>  $ estimator.type  : chr "classical"
#>  $ n.data          : int 30
#>  $ lambda          : num 1
#>  $ trend           : chr "cte"
#>  $ pairs.min       : num 5
#>  $ nugget.tolerance: num 1e-12
#>  $ direction       : chr "omnidirectional"
#>  $ tolerance       : chr "none"
#>  $ uvec            : num [1:30] 1 2.19 2.61 3.11 3.7 ...
#>  $ call            : language variog(geodata = geodata, breaks = breaks, trend = trend, pairs.min = 5)
#> 
#> vm
#> 'variomodel', variofit'
#> List of 17
#>  $ nugget               : num 0.931
#>  $ cov.pars             : num [1:2] 0 160
#>  $ cov.model            : chr "gaussian"
#>  $ kappa                : num 0.5
#>  $ value                : num 19.9
#>  $ trend                : chr "cte"
#>  $ beta.ols             : num -5.08e-09
#>  $ practicalRange       : num 277
#>  $ max.dist             : num 320
#>  $ minimisation.function: chr "optim"
#>  $ weights              : chr "npairs"
#>  $ method               : chr "WLS"
#>  $ fix.nugget           : logi FALSE
#>  $ fix.kappa            : logi TRUE
#>  $ lambda               : num 1
#>  $ message              : chr "optim convergence code: 0"
#>  $ call                 : language variofit(vario = vg, ini.cov.pars = c(initialVal, startRange), cov.model = varModels[i],      nugget = initialVal)
#> 
#> var: c 
#> df
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1150 obs. of  3 variables:
#>  $ x: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ y: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z: num  0.837 0.811 0.785 0.759 0.733 ...
#> 
#> df.poly
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1150 obs. of  3 variables:
#>  $ gx: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ gy: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z : num  0.837 0.811 0.785 0.759 0.733 ...
#> 
#> lambda
#> 'numeric'
#>  num 1
#> 
#> vg
#> 'variogram'
#> List of 20
#>  $ u               : num [1:9] 60.9 86.5 103 122.7 146.1 ...
#>  $ v               : num [1:9] 0.834 1.239 2.593 1.597 0.621 ...
#>  $ n               : num [1:9] 7 9 10 10 18 19 36 34 38
#>  $ sd              : num [1:9] 1.22 1.41 1.86 1.47 1.19 ...
#>  $ bins.lim        : num [1:31] 1.00e-12 2.00 2.38 2.84 3.38 ...
#>  $ ind.bin         : logi [1:30] FALSE FALSE FALSE FALSE FALSE FALSE ...
#>  $ var.mark        : num 0.931
#>  $ beta.ols        : num 1.77e-11
#>  $ output.type     : chr "bin"
#>  $ max.dist        : num 320
#>  $ estimator.type  : chr "classical"
#>  $ n.data          : int 30
#>  $ lambda          : num 1
#>  $ trend           : chr "cte"
#>  $ pairs.min       : num 5
#>  $ nugget.tolerance: num 1e-12
#>  $ direction       : chr "omnidirectional"
#>  $ tolerance       : chr "none"
#>  $ uvec            : num [1:30] 1 2.19 2.61 3.11 3.7 ...
#>  $ call            : language variog(geodata = geodata, breaks = breaks, trend = trend, pairs.min = 5)
#> 
#> vm
#> 'variomodel', variofit'
#> List of 17
#>  $ nugget               : num 1.04
#>  $ cov.pars             : num [1:2] 0 161
#>  $ cov.model            : chr "exponential"
#>  $ kappa                : num 0.5
#>  $ value                : num 40.1
#>  $ trend                : chr "cte"
#>  $ beta.ols             : num 1.77e-11
#>  $ practicalRange       : num 481
#>  $ max.dist             : num 320
#>  $ minimisation.function: chr "optim"
#>  $ weights              : chr "npairs"
#>  $ method               : chr "WLS"
#>  $ fix.nugget           : logi FALSE
#>  $ fix.kappa            : logi TRUE
#>  $ lambda               : num 1
#>  $ message              : chr "optim convergence code: 0"
#>  $ call                 : language variofit(vario = vg, ini.cov.pars = c(initialVal, startRange), cov.model = varModels[i],      nugget = initialVal)
#> 
#> var: p 
#> df
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1150 obs. of  3 variables:
#>  $ x: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ y: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z: num  0.375 0.337 0.301 0.266 0.233 ...
#> 
#> df.poly
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1150 obs. of  3 variables:
#>  $ gx: num  10 30 50 70 90 110 130 150 170 190 ...
#>  $ gy: num  10 10 10 10 10 10 10 10 10 10 ...
#>  $ z : num  0.375 0.337 0.301 0.266 0.233 ...
#> 
#> lambda
#> 'numeric'
#>  num 1
#> 
#> vg
#> 'variogram'
#> List of 20
#>  $ u               : num [1:9] 60.9 86.5 103 122.7 146.1 ...
#>  $ v               : num [1:9] 1.311 1.325 0.366 1.332 1.275 ...
#>  $ n               : num [1:9] 7 9 10 10 18 19 36 34 38
#>  $ sd              : num [1:9] 1.963 1.371 0.441 1.355 1.616 ...
#>  $ bins.lim        : num [1:31] 1.00e-12 2.00 2.38 2.84 3.38 ...
#>  $ ind.bin         : logi [1:30] FALSE FALSE FALSE FALSE FALSE FALSE ...
#>  $ var.mark        : num 0.886
#>  $ beta.ols        : num -1.19e-09
#>  $ output.type     : chr "bin"
#>  $ max.dist        : num 320
#>  $ estimator.type  : chr "classical"
#>  $ n.data          : int 30
#>  $ lambda          : num 1
#>  $ trend           : chr "cte"
#>  $ pairs.min       : num 5
#>  $ nugget.tolerance: num 1e-12
#>  $ direction       : chr "omnidirectional"
#>  $ tolerance       : chr "none"
#>  $ uvec            : num [1:30] 1 2.19 2.61 3.11 3.7 ...
#>  $ call            : language variog(geodata = geodata, breaks = breaks, trend = trend, pairs.min = 5)
#> 
#> vm
#> 'variomodel', variofit'
#> List of 17
#>  $ nugget               : num 1.01
#>  $ cov.pars             : num [1:2] 0 160
#>  $ cov.model            : chr "circular"
#>  $ kappa                : num 0.5
#>  $ value                : num 8.99
#>  $ trend                : chr "cte"
#>  $ beta.ols             : num -1.19e-09
#>  $ practicalRange       : num 160
#>  $ max.dist             : num 320
#>  $ minimisation.function: chr "optim"
#>  $ weights              : chr "npairs"
#>  $ method               : chr "WLS"
#>  $ fix.nugget           : logi FALSE
#>  $ fix.kappa            : logi TRUE
#>  $ lambda               : num 1
#>  $ message              : chr "optim convergence code: 0"
#>  $ call                 : language variofit(vario = vg, ini.cov.pars = c(initialVal, startRange), cov.model = varModels[i],      nugget = initialVal)

Transform krig results to a tibble for easier manipulation and visualization.

Visualize multiple soil variables

The structure of this dataframe allows faceting with gglot2.

ggplot(krig_df, aes(x, y, fill = z)) +
  geom_tile() +
  facet_wrap("var") +
  coord_equal()