R Under development (unstable) (2025-02-12 r87715 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > options(digits=7) > > library(automap) > library(sp) > # Neccessary to silence sf startup messages > suppressMessages(library(sf)) > > > # Data preparation > data(meuse) > coordinates(meuse) =~ x+y > data(meuse.grid) > gridded(meuse.grid) =~ x+y > > # Fitting some variograms > variogram = autofitVariogram(zinc ~ soil + ffreq + dist, meuse, + miscFitOptions = list(min.np.bin = 500)) > variogram $exp_var np dist gamma dir.hor dir.ver id 1 500 174.0355 38065.17 0 0 var1 2 711 337.6036 50150.73 0 0 var1 3 830 502.0477 50347.61 0 0 var1 4 1349 713.2149 61255.29 0 0 var1 5 1314 961.2718 70311.91 0 0 var1 6 1139 1213.4116 75870.33 0 0 var1 7 1355 1506.5505 72043.88 0 0 var1 $var_model model psill range kappa 1 Nug 5294.22 0.000 0.0 2 Ste 91418.10 1781.645 0.2 $sserr [1] 130690.3 attr(,"class") [1] "autofitVariogram" "list" > > # ...and diable the merging, note the difference between the two plots > variogram = autofitVariogram(zinc ~ soil + ffreq + dist, meuse, + miscFitOptions = list(min.np.bin = 500, merge.small.bins = FALSE)) > variogram $exp_var np dist gamma dir.hor dir.ver id 1 17 59.33470 23959.65 0 0 var1 2 36 86.01449 16671.65 0 0 var1 3 114 131.02870 33270.31 0 0 var1 4 149 176.18845 52846.33 0 0 var1 5 184 226.75652 34555.29 0 0 var1 6 711 337.60359 50150.73 0 0 var1 7 830 502.04773 50347.61 0 0 var1 8 1349 713.21485 61255.29 0 0 var1 9 1314 961.27179 70311.91 0 0 var1 10 1139 1213.41157 75870.33 0 0 var1 11 1355 1506.55052 72043.88 0 0 var1 $var_model model psill range kappa 1 Nug 2454.939 0.0000 0.0 2 Ste 68746.147 482.0245 0.3 $sserr [1] 2044182 attr(,"class") [1] "autofitVariogram" "list" > > # Ordinary kriging > kriging_result = autoKrige(zinc~1, meuse, meuse.grid) [using ordinary kriging] > summary(kriging_result) krige_output: Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: var1.pred var1.var var1.stdev Min. : 106.1 Min. : 40440 Min. :201.1 1st Qu.: 207.1 1st Qu.: 49617 1st Qu.:222.7 Median : 309.4 Median : 55256 Median :235.1 Mean : 406.2 Mean : 62453 Mean :247.5 3rd Qu.: 555.1 3rd Qu.: 69014 3rd Qu.:262.7 Max. :1620.1 Max. :142357 Max. :377.3 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 46207.03 0 0 var1 2 36 86.01449 32203.57 0 0 var1 3 114 131.02870 53862.70 0 0 var1 4 149 176.18845 87068.48 0 0 var1 5 184 226.75652 62486.31 0 0 var1 6 711 337.60359 104453.77 0 0 var1 7 830 502.04773 127330.49 0 0 var1 8 1349 713.21485 150556.52 0 0 var1 9 1314 961.27179 165153.58 0 0 var1 10 1139 1213.41157 165495.97 0 0 var1 11 1355 1506.55052 149245.86 0 0 var1 var_model: model psill range kappa 1 Nug 31957.6 0.0000 0.0 2 Ste 133556.2 485.5114 1.1 Sums of squares betw. var. model and sample var.[1] 4235991 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value:Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: upper lower position Min. : 554.1 Min. :-448.11 higher : 717 1st Qu.: 673.3 1st Qu.:-272.23 not distinguishable:2386 Median : 812.2 Median :-156.60 Mean : 891.3 Mean : -78.90 3rd Qu.:1051.5 3rd Qu.: 50.58 Max. :2065.8 Max. :1174.50 $p [1] 95 $value [1] 75 > > kriging_result = autoKrige(zinc~1, meuse, meuse.grid, fix.values = c(0.2,NA,NA)) [using ordinary kriging] > summary(kriging_result) krige_output: Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: var1.pred var1.var var1.stdev Min. : 119.8 Min. : 3295 Min. : 57.4 1st Qu.: 201.3 1st Qu.: 40754 1st Qu.:201.9 Median : 298.4 Median : 52570 Median :229.3 Mean : 403.7 Mean : 58013 Mean :235.4 3rd Qu.: 552.0 3rd Qu.: 71373 3rd Qu.:267.2 Max. :1800.9 Max. :141490 Max. :376.2 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 46207.03 0 0 var1 2 36 86.01449 32203.57 0 0 var1 3 114 131.02870 53862.70 0 0 var1 4 149 176.18845 87068.48 0 0 var1 5 184 226.75652 62486.31 0 0 var1 6 711 337.60359 104453.77 0 0 var1 7 830 502.04773 127330.49 0 0 var1 8 1349 713.21485 150556.52 0 0 var1 9 1314 961.27179 165153.58 0 0 var1 10 1139 1213.41157 165495.97 0 0 var1 11 1355 1506.55052 149245.86 0 0 var1 var_model: model psill range kappa 1 Nug 0.2 0.0000 0.0 2 Ste 173889.1 567.4229 0.4 Sums of squares betw. var. model and sample var.[1] 4751276 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value:Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: upper lower position Min. : 292.2 Min. :-465.4 higher : 759 1st Qu.: 654.8 1st Qu.:-260.3 not distinguishable:2344 Median : 788.4 Median :-130.9 Mean : 865.1 Mean : -57.8 3rd Qu.:1027.8 3rd Qu.: 70.0 Max. :2108.8 Max. :1575.5 $p [1] 95 $value [1] 75 > > # Universal kriging > kriging_result = autoKrige(zinc~soil+ffreq+dist, meuse, meuse.grid) [using universal kriging] > summary(kriging_result) krige_output: Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: var1.pred var1.var var1.stdev Min. :-181.0 Min. : 9538 Min. : 97.66 1st Qu.: 114.4 1st Qu.:29049 1st Qu.:170.44 Median : 258.9 Median :34707 Median :186.30 Mean : 344.9 Mean :36358 Mean :188.75 3rd Qu.: 520.9 3rd Qu.:42945 3rd Qu.:207.23 Max. :1727.5 Max. :67433 Max. :259.68 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 23959.65 0 0 var1 2 36 86.01449 16671.65 0 0 var1 3 114 131.02870 33270.31 0 0 var1 4 149 176.18845 52846.33 0 0 var1 5 184 226.75652 34555.29 0 0 var1 6 711 337.60359 50150.73 0 0 var1 7 830 502.04773 50347.61 0 0 var1 8 1349 713.21485 61255.29 0 0 var1 9 1314 961.27179 70311.91 0 0 var1 10 1139 1213.41157 75870.33 0 0 var1 11 1355 1506.55052 72043.88 0 0 var1 var_model: model psill range kappa 1 Nug 2454.939 0.0000 0.0 2 Ste 68746.147 482.0245 0.3 Sums of squares betw. var. model and sample var.[1] 2044182 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value:Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: upper lower position Min. : 166.2 Min. :-662.81 higher : 933 1st Qu.: 481.0 1st Qu.:-248.96 not distinguishable:2170 Median : 609.8 Median : -90.68 Mean : 714.8 Mean : -25.09 3rd Qu.: 912.3 3rd Qu.: 140.09 Max. :1969.6 Max. :1485.50 $p [1] 95 $value [1] 75 > > # Block kriging > kriging_result = autoKrige(zinc~soil+ffreq+dist, meuse, meuse.grid, block = c(400,400)) [using universal kriging] > summary(kriging_result) krige_output: Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: var1.pred var1.var var1.stdev Min. :-190.2 Min. : 2251 Min. : 47.45 1st Qu.: 111.5 1st Qu.: 4680 1st Qu.: 68.41 Median : 283.5 Median : 6688 Median : 81.78 Mean : 348.1 Mean : 8516 Mean : 88.65 3rd Qu.: 535.6 3rd Qu.:10855 3rd Qu.:104.19 Max. :1338.2 Max. :29401 Max. :171.47 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 23959.65 0 0 var1 2 36 86.01449 16671.65 0 0 var1 3 114 131.02870 33270.31 0 0 var1 4 149 176.18845 52846.33 0 0 var1 5 184 226.75652 34555.29 0 0 var1 6 711 337.60359 50150.73 0 0 var1 7 830 502.04773 50347.61 0 0 var1 8 1349 713.21485 61255.29 0 0 var1 9 1314 961.27179 70311.91 0 0 var1 10 1139 1213.41157 75870.33 0 0 var1 11 1355 1506.55052 72043.88 0 0 var1 var_model: model psill range kappa 1 Nug 2454.939 0.0000 0.0 2 Ste 68746.147 482.0245 0.3 Sums of squares betw. var. model and sample var.[1] 2044182 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value:Object of class SpatialPixelsDataFrame Coordinates: min max x 178440 181560 y 329600 333760 Is projected: NA proj4string : [NA] Number of points: 3103 Grid attributes: cellcentre.offset cellsize cells.dim x 178460 40 78 y 329620 40 104 Data attributes: upper lower position Min. : 23.99 Min. :-490.57 higher :1790 1st Qu.: 276.97 1st Qu.: -47.45 lower : 107 Median : 433.15 Median : 131.75 not distinguishable:1206 Mean : 521.89 Mean : 174.40 3rd Qu.: 732.91 3rd Qu.: 359.36 Max. :1480.03 Max. :1196.43 $p [1] 95 $value [1] 75 > > # Kriging with power model > kriging_result = autoKrige(zinc~soil+ffreq+dist, meuse, meuse.grid, model = "Pow") [using universal kriging] Warning message: In getModel(initial_sill - initial_nugget, m, initial_range, kappa = 0, : Using the power model is at your own risk, read the docs of autofitVariogram for more details. > > > > # Testing with sf objects > meuse = as(meuse, "sf") > meuse.grid = as(meuse.grid, "sf") > > # Fitting some variograms > variogram = autofitVariogram(zinc ~ soil + ffreq + dist, meuse, + miscFitOptions = list(min.np.bin = 500)) > variogram $exp_var np dist gamma dir.hor dir.ver id 1 500 174.0355 38065.17 0 0 var1 2 711 337.6036 50150.73 0 0 var1 3 830 502.0477 50347.61 0 0 var1 4 1349 713.2149 61255.29 0 0 var1 5 1314 961.2718 70311.91 0 0 var1 6 1139 1213.4116 75870.33 0 0 var1 7 1355 1506.5505 72043.88 0 0 var1 $var_model model psill range kappa 1 Nug 5294.22 0.000 0.0 2 Ste 91418.10 1781.645 0.2 $sserr [1] 130690.3 attr(,"class") [1] "autofitVariogram" "list" > > # ...and diable the merging, note the difference between the two plots > variogram = autofitVariogram(zinc ~ soil + ffreq + dist, meuse, + miscFitOptions = list(min.np.bin = 500, merge.small.bins = FALSE)) > variogram $exp_var np dist gamma dir.hor dir.ver id 1 17 59.33470 23959.65 0 0 var1 2 36 86.01449 16671.65 0 0 var1 3 114 131.02870 33270.31 0 0 var1 4 149 176.18845 52846.33 0 0 var1 5 184 226.75652 34555.29 0 0 var1 6 711 337.60359 50150.73 0 0 var1 7 830 502.04773 50347.61 0 0 var1 8 1349 713.21485 61255.29 0 0 var1 9 1314 961.27179 70311.91 0 0 var1 10 1139 1213.41157 75870.33 0 0 var1 11 1355 1506.55052 72043.88 0 0 var1 $var_model model psill range kappa 1 Nug 2454.939 0.0000 0.0 2 Ste 68746.147 482.0245 0.3 $sserr [1] 2044182 attr(,"class") [1] "autofitVariogram" "list" > > # Ordinary kriging > kriging_result = autoKrige(zinc~1, meuse, meuse.grid) [using ordinary kriging] > summary(kriging_result) krige_output: var1.pred var1.var geometry var1.stdev Min. : 106.1 Min. : 40440 POINT :3103 Min. :201.1 1st Qu.: 207.1 1st Qu.: 49617 epsg:NA: 0 1st Qu.:222.7 Median : 309.4 Median : 55256 Median :235.1 Mean : 406.2 Mean : 62453 Mean :247.5 3rd Qu.: 555.1 3rd Qu.: 69014 3rd Qu.:262.7 Max. :1620.1 Max. :142357 Max. :377.3 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 46207.03 0 0 var1 2 36 86.01449 32203.57 0 0 var1 3 114 131.02870 53862.70 0 0 var1 4 149 176.18845 87068.48 0 0 var1 5 184 226.75652 62486.31 0 0 var1 6 711 337.60359 104453.77 0 0 var1 7 830 502.04773 127330.49 0 0 var1 8 1349 713.21485 150556.52 0 0 var1 9 1314 961.27179 165153.58 0 0 var1 10 1139 1213.41157 165495.97 0 0 var1 11 1355 1506.55052 149245.86 0 0 var1 var_model: model psill range kappa 1 Nug 31957.6 0.0000 0.0 2 Ste 133556.2 485.5114 1.1 Sums of squares betw. var. model and sample var.[1] 4235991 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position geometry Min. : 554.1 Min. :-448.11 higher : 717 POINT :3103 1st Qu.: 673.3 1st Qu.:-272.23 not distinguishable:2386 epsg:NA: 0 Median : 812.2 Median :-156.60 Mean : 891.3 Mean : -78.90 3rd Qu.:1051.5 3rd Qu.: 50.58 Max. :2065.8 Max. :1174.50 $p [1] 95 $value [1] 75 > > kriging_result = autoKrige(zinc~1, meuse, meuse.grid, fix.values = c(0.2,NA,NA)) [using ordinary kriging] > summary(kriging_result) krige_output: var1.pred var1.var geometry var1.stdev Min. : 119.8 Min. : 3295 POINT :3103 Min. : 57.4 1st Qu.: 201.3 1st Qu.: 40754 epsg:NA: 0 1st Qu.:201.9 Median : 298.4 Median : 52570 Median :229.3 Mean : 403.7 Mean : 58013 Mean :235.4 3rd Qu.: 552.0 3rd Qu.: 71373 3rd Qu.:267.2 Max. :1800.9 Max. :141490 Max. :376.2 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 46207.03 0 0 var1 2 36 86.01449 32203.57 0 0 var1 3 114 131.02870 53862.70 0 0 var1 4 149 176.18845 87068.48 0 0 var1 5 184 226.75652 62486.31 0 0 var1 6 711 337.60359 104453.77 0 0 var1 7 830 502.04773 127330.49 0 0 var1 8 1349 713.21485 150556.52 0 0 var1 9 1314 961.27179 165153.58 0 0 var1 10 1139 1213.41157 165495.97 0 0 var1 11 1355 1506.55052 149245.86 0 0 var1 var_model: model psill range kappa 1 Nug 0.2 0.0000 0.0 2 Ste 173889.1 567.4229 0.4 Sums of squares betw. var. model and sample var.[1] 4751276 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position geometry Min. : 292.2 Min. :-465.4 higher : 759 POINT :3103 1st Qu.: 654.8 1st Qu.:-260.3 not distinguishable:2344 epsg:NA: 0 Median : 788.4 Median :-130.9 Mean : 865.1 Mean : -57.8 3rd Qu.:1027.8 3rd Qu.: 70.0 Max. :2108.8 Max. :1575.5 $p [1] 95 $value [1] 75 > > # Universal kriging > kriging_result = autoKrige(zinc~soil+ffreq+dist, meuse, meuse.grid) [using universal kriging] > summary(kriging_result) krige_output: var1.pred var1.var geometry var1.stdev Min. :-181.0 Min. : 9538 POINT :3103 Min. : 97.66 1st Qu.: 114.4 1st Qu.:29049 epsg:NA: 0 1st Qu.:170.44 Median : 258.9 Median :34707 Median :186.30 Mean : 344.9 Mean :36358 Mean :188.75 3rd Qu.: 520.9 3rd Qu.:42945 3rd Qu.:207.23 Max. :1727.5 Max. :67433 Max. :259.68 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 23959.65 0 0 var1 2 36 86.01449 16671.65 0 0 var1 3 114 131.02870 33270.31 0 0 var1 4 149 176.18845 52846.33 0 0 var1 5 184 226.75652 34555.29 0 0 var1 6 711 337.60359 50150.73 0 0 var1 7 830 502.04773 50347.61 0 0 var1 8 1349 713.21485 61255.29 0 0 var1 9 1314 961.27179 70311.91 0 0 var1 10 1139 1213.41157 75870.33 0 0 var1 11 1355 1506.55052 72043.88 0 0 var1 var_model: model psill range kappa 1 Nug 2454.939 0.0000 0.0 2 Ste 68746.147 482.0245 0.3 Sums of squares betw. var. model and sample var.[1] 2044182 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position geometry Min. : 166.2 Min. :-662.81 higher : 933 POINT :3103 1st Qu.: 481.0 1st Qu.:-248.96 not distinguishable:2170 epsg:NA: 0 Median : 609.8 Median : -90.68 Mean : 714.8 Mean : -25.09 3rd Qu.: 912.3 3rd Qu.: 140.09 Max. :1969.6 Max. :1485.50 $p [1] 95 $value [1] 75 > > # Block kriging > kriging_result = autoKrige(zinc~soil+ffreq+dist, meuse, meuse.grid, block = c(400,400)) [using universal kriging] > summary(kriging_result) krige_output: var1.pred var1.var geometry var1.stdev Min. :-190.2 Min. : 2251 POINT :3103 Min. : 47.45 1st Qu.: 111.5 1st Qu.: 4680 epsg:NA: 0 1st Qu.: 68.41 Median : 283.5 Median : 6688 Median : 81.78 Mean : 348.1 Mean : 8516 Mean : 88.65 3rd Qu.: 535.6 3rd Qu.:10855 3rd Qu.:104.19 Max. :1338.2 Max. :29401 Max. :171.47 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.33470 23959.65 0 0 var1 2 36 86.01449 16671.65 0 0 var1 3 114 131.02870 33270.31 0 0 var1 4 149 176.18845 52846.33 0 0 var1 5 184 226.75652 34555.29 0 0 var1 6 711 337.60359 50150.73 0 0 var1 7 830 502.04773 50347.61 0 0 var1 8 1349 713.21485 61255.29 0 0 var1 9 1314 961.27179 70311.91 0 0 var1 10 1139 1213.41157 75870.33 0 0 var1 11 1355 1506.55052 72043.88 0 0 var1 var_model: model psill range kappa 1 Nug 2454.939 0.0000 0.0 2 Ste 68746.147 482.0245 0.3 Sums of squares betw. var. model and sample var.[1] 2044182 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position geometry Min. : 23.99 Min. :-490.57 higher :1790 POINT :3103 1st Qu.: 276.97 1st Qu.: -47.45 lower : 107 epsg:NA: 0 Median : 433.15 Median : 131.75 not distinguishable:1206 Mean : 521.89 Mean : 174.40 3rd Qu.: 732.91 3rd Qu.: 359.36 Max. :1480.03 Max. :1196.43 $p [1] 95 $value [1] 75 > > # Kriging with power model > kriging_result = autoKrige(zinc~soil+ffreq+dist, meuse, meuse.grid, model = "Pow") [using universal kriging] Warning message: In getModel(initial_sill - initial_nugget, m, initial_range, kappa = 0, : Using the power model is at your own risk, read the docs of autofitVariogram for more details. > > > proc.time() user system elapsed 15.21 0.34 15.54