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=8) > > # Silence version differences > suppressMessages(library(automap)) > suppressMessages(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.03552 38065.172 0 0 var1 2 711 337.60359 50150.725 0 0 var1 3 830 502.04773 50347.613 0 0 var1 4 1349 713.21485 61255.291 0 0 var1 5 1314 961.27179 70311.908 0 0 var1 6 1139 1213.41157 75870.335 0 0 var1 7 1355 1506.55052 72043.880 0 0 var1 $var_model model psill range kappa 1 Nug 5294.2204 0.0000 0.0 2 Ste 91418.0987 1781.6454 0.2 $sserr [1] 130690.32 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.334696 23959.646 0 0 var1 2 36 86.014494 16671.647 0 0 var1 3 114 131.028697 33270.313 0 0 var1 4 149 176.188454 52846.333 0 0 var1 5 184 226.756517 34555.290 0 0 var1 6 711 337.603595 50150.725 0 0 var1 7 830 502.047734 50347.613 0 0 var1 8 1349 713.214854 61255.291 0 0 var1 9 1314 961.271787 70311.908 0 0 var1 10 1139 1213.411568 75870.335 0 0 var1 11 1355 1506.550518 72043.880 0 0 var1 $var_model model psill range kappa 1 Nug 2454.9391 0.00000 0.0 2 Ste 68746.1473 482.02453 0.3 $sserr [1] 2044181.5 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.06 Min. : 40440 Min. :201.10 1st Qu.: 207.08 1st Qu.: 49617 1st Qu.:222.75 Median : 309.35 Median : 55256 Median :235.07 Mean : 406.18 Mean : 62453 Mean :247.50 3rd Qu.: 555.12 3rd Qu.: 69014 3rd Qu.:262.70 Max. :1620.15 Max. :142357 Max. :377.30 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 46207.029 0 0 var1 2 36 86.014494 32203.569 0 0 var1 3 114 131.028697 53862.702 0 0 var1 4 149 176.188454 87068.477 0 0 var1 5 184 226.756517 62486.307 0 0 var1 6 711 337.603595 104453.774 0 0 var1 7 830 502.047734 127330.495 0 0 var1 8 1349 713.214854 150556.520 0 0 var1 9 1314 961.271787 165153.577 0 0 var1 10 1139 1213.411568 165495.972 0 0 var1 11 1355 1506.550518 149245.856 0 0 var1 var_model: model psill range kappa 1 Nug 31957.60 0.00000 0.0 2 Ste 133556.19 485.51142 1.1 Sums of squares betw. var. model and sample var.[1] 4235990.6 > 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.06 Min. :-448.109 higher : 717 1st Qu.: 673.32 1st Qu.:-272.229 not distinguishable:2386 Median : 812.24 Median :-156.599 Mean : 891.27 Mean : -78.905 3rd Qu.:1051.51 3rd Qu.: 50.584 Max. :2065.79 Max. :1174.502 $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.83 Min. : 3294.8 Min. : 57.40 1st Qu.: 201.35 1st Qu.: 40753.6 1st Qu.:201.88 Median : 298.36 Median : 52570.2 Median :229.28 Mean : 403.67 Mean : 58013.3 Mean :235.45 3rd Qu.: 551.98 3rd Qu.: 71373.5 3rd Qu.:267.16 Max. :1800.86 Max. :141489.6 Max. :376.15 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 46207.029 0 0 var1 2 36 86.014494 32203.569 0 0 var1 3 114 131.028697 53862.702 0 0 var1 4 149 176.188454 87068.477 0 0 var1 5 184 226.756517 62486.307 0 0 var1 6 711 337.603595 104453.774 0 0 var1 7 830 502.047734 127330.495 0 0 var1 8 1349 713.214854 150556.520 0 0 var1 9 1314 961.271787 165153.577 0 0 var1 10 1139 1213.411568 165495.972 0 0 var1 11 1355 1506.550518 149245.856 0 0 var1 var_model: model psill range kappa 1 Nug 0.20 0.00000 0.0 2 Ste 173889.06 567.42291 0.4 Sums of squares betw. var. model and sample var.[1] 4751276.5 > 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.16 Min. :-465.396 higher : 759 1st Qu.: 654.80 1st Qu.:-260.298 not distinguishable:2344 Median : 788.39 Median :-130.853 Mean : 865.14 Mean : -57.804 3rd Qu.:1027.82 3rd Qu.: 70.000 Max. :2108.76 Max. :1575.547 $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. :-180.97 Min. : 9538 Min. : 97.663 1st Qu.: 114.37 1st Qu.:29049 1st Qu.:170.439 Median : 258.94 Median :34707 Median :186.298 Mean : 344.85 Mean :36358 Mean :188.749 3rd Qu.: 520.93 3rd Qu.:42945 3rd Qu.:207.231 Max. :1727.54 Max. :67433 Max. :259.680 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 23959.646 0 0 var1 2 36 86.014494 16671.647 0 0 var1 3 114 131.028697 33270.313 0 0 var1 4 149 176.188454 52846.333 0 0 var1 5 184 226.756517 34555.290 0 0 var1 6 711 337.603595 50150.725 0 0 var1 7 830 502.047734 50347.613 0 0 var1 8 1349 713.214854 61255.291 0 0 var1 9 1314 961.271787 70311.908 0 0 var1 10 1139 1213.411568 75870.335 0 0 var1 11 1355 1506.550518 72043.880 0 0 var1 var_model: model psill range kappa 1 Nug 2454.9391 0.00000 0.0 2 Ste 68746.1473 482.02453 0.3 Sums of squares betw. var. model and sample var.[1] 2044181.5 > 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.19 Min. :-662.805 higher : 933 1st Qu.: 480.97 1st Qu.:-248.957 not distinguishable:2170 Median : 609.77 Median : -90.682 Mean : 714.79 Mean : -25.087 3rd Qu.: 912.27 3rd Qu.: 140.087 Max. :1969.57 Max. :1485.503 $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.15 Min. : 2251.5 Min. : 47.450 1st Qu.: 111.52 1st Qu.: 4679.9 1st Qu.: 68.410 Median : 283.45 Median : 6688.4 Median : 81.783 Mean : 348.14 Mean : 8516.3 Mean : 88.646 3rd Qu.: 535.64 3rd Qu.:10854.7 3rd Qu.:104.186 Max. :1338.23 Max. :29401.2 Max. :171.468 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 23959.646 0 0 var1 2 36 86.014494 16671.647 0 0 var1 3 114 131.028697 33270.313 0 0 var1 4 149 176.188454 52846.333 0 0 var1 5 184 226.756517 34555.290 0 0 var1 6 711 337.603595 50150.725 0 0 var1 7 830 502.047734 50347.613 0 0 var1 8 1349 713.214854 61255.291 0 0 var1 9 1314 961.271787 70311.908 0 0 var1 10 1139 1213.411568 75870.335 0 0 var1 11 1355 1506.550518 72043.880 0 0 var1 var_model: model psill range kappa 1 Nug 2454.9391 0.00000 0.0 2 Ste 68746.1473 482.02453 0.3 Sums of squares betw. var. model and sample var.[1] 2044181.5 > 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.995 Min. :-490.572 higher :1790 1st Qu.: 276.975 1st Qu.: -47.453 lower : 107 Median : 433.154 Median : 131.751 not distinguishable:1206 Mean : 521.886 Mean : 174.399 3rd Qu.: 732.911 3rd Qu.: 359.361 Max. :1480.025 Max. :1196.429 $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.03552 38065.172 0 0 var1 2 711 337.60359 50150.725 0 0 var1 3 830 502.04773 50347.613 0 0 var1 4 1349 713.21485 61255.291 0 0 var1 5 1314 961.27179 70311.908 0 0 var1 6 1139 1213.41157 75870.335 0 0 var1 7 1355 1506.55052 72043.880 0 0 var1 $var_model model psill range kappa 1 Nug 5294.2204 0.0000 0.0 2 Ste 91418.0987 1781.6454 0.2 $sserr [1] 130690.32 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.334696 23959.646 0 0 var1 2 36 86.014494 16671.647 0 0 var1 3 114 131.028697 33270.313 0 0 var1 4 149 176.188454 52846.333 0 0 var1 5 184 226.756517 34555.290 0 0 var1 6 711 337.603595 50150.725 0 0 var1 7 830 502.047734 50347.613 0 0 var1 8 1349 713.214854 61255.291 0 0 var1 9 1314 961.271787 70311.908 0 0 var1 10 1139 1213.411568 75870.335 0 0 var1 11 1355 1506.550518 72043.880 0 0 var1 $var_model model psill range kappa 1 Nug 2454.9391 0.00000 0.0 2 Ste 68746.1473 482.02453 0.3 $sserr [1] 2044181.5 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.06 Min. : 40440 POINT :3103 Min. :201.10 1st Qu.: 207.08 1st Qu.: 49617 epsg:NA: 0 1st Qu.:222.75 Median : 309.35 Median : 55256 Median :235.07 Mean : 406.18 Mean : 62453 Mean :247.50 3rd Qu.: 555.12 3rd Qu.: 69014 3rd Qu.:262.70 Max. :1620.15 Max. :142357 Max. :377.30 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 46207.029 0 0 var1 2 36 86.014494 32203.569 0 0 var1 3 114 131.028697 53862.702 0 0 var1 4 149 176.188454 87068.477 0 0 var1 5 184 226.756517 62486.307 0 0 var1 6 711 337.603595 104453.774 0 0 var1 7 830 502.047734 127330.495 0 0 var1 8 1349 713.214854 150556.520 0 0 var1 9 1314 961.271787 165153.577 0 0 var1 10 1139 1213.411568 165495.972 0 0 var1 11 1355 1506.550518 149245.856 0 0 var1 var_model: model psill range kappa 1 Nug 31957.60 0.00000 0.0 2 Ste 133556.19 485.51142 1.1 Sums of squares betw. var. model and sample var.[1] 4235990.6 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position geometry Min. : 554.06 Min. :-448.109 higher : 717 POINT :3103 1st Qu.: 673.32 1st Qu.:-272.229 not distinguishable:2386 epsg:NA: 0 Median : 812.24 Median :-156.599 Mean : 891.27 Mean : -78.905 3rd Qu.:1051.51 3rd Qu.: 50.584 Max. :2065.79 Max. :1174.502 $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.83 Min. : 3294.8 POINT :3103 Min. : 57.40 1st Qu.: 201.35 1st Qu.: 40753.6 epsg:NA: 0 1st Qu.:201.88 Median : 298.36 Median : 52570.2 Median :229.28 Mean : 403.67 Mean : 58013.3 Mean :235.45 3rd Qu.: 551.98 3rd Qu.: 71373.5 3rd Qu.:267.16 Max. :1800.86 Max. :141489.6 Max. :376.15 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 46207.029 0 0 var1 2 36 86.014494 32203.569 0 0 var1 3 114 131.028697 53862.702 0 0 var1 4 149 176.188454 87068.477 0 0 var1 5 184 226.756517 62486.307 0 0 var1 6 711 337.603595 104453.774 0 0 var1 7 830 502.047734 127330.495 0 0 var1 8 1349 713.214854 150556.520 0 0 var1 9 1314 961.271787 165153.577 0 0 var1 10 1139 1213.411568 165495.972 0 0 var1 11 1355 1506.550518 149245.856 0 0 var1 var_model: model psill range kappa 1 Nug 0.20 0.00000 0.0 2 Ste 173889.06 567.42291 0.4 Sums of squares betw. var. model and sample var.[1] 4751276.5 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position geometry Min. : 292.16 Min. :-465.396 higher : 759 POINT :3103 1st Qu.: 654.80 1st Qu.:-260.298 not distinguishable:2344 epsg:NA: 0 Median : 788.39 Median :-130.853 Mean : 865.14 Mean : -57.804 3rd Qu.:1027.82 3rd Qu.: 70.000 Max. :2108.76 Max. :1575.547 $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. :-180.97 Min. : 9538 POINT :3103 Min. : 97.663 1st Qu.: 114.37 1st Qu.:29049 epsg:NA: 0 1st Qu.:170.439 Median : 258.94 Median :34707 Median :186.298 Mean : 344.85 Mean :36358 Mean :188.749 3rd Qu.: 520.93 3rd Qu.:42945 3rd Qu.:207.231 Max. :1727.54 Max. :67433 Max. :259.680 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 23959.646 0 0 var1 2 36 86.014494 16671.647 0 0 var1 3 114 131.028697 33270.313 0 0 var1 4 149 176.188454 52846.333 0 0 var1 5 184 226.756517 34555.290 0 0 var1 6 711 337.603595 50150.725 0 0 var1 7 830 502.047734 50347.613 0 0 var1 8 1349 713.214854 61255.291 0 0 var1 9 1314 961.271787 70311.908 0 0 var1 10 1139 1213.411568 75870.335 0 0 var1 11 1355 1506.550518 72043.880 0 0 var1 var_model: model psill range kappa 1 Nug 2454.9391 0.00000 0.0 2 Ste 68746.1473 482.02453 0.3 Sums of squares betw. var. model and sample var.[1] 2044181.5 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position geometry Min. : 166.19 Min. :-662.805 higher : 933 POINT :3103 1st Qu.: 480.97 1st Qu.:-248.957 not distinguishable:2170 epsg:NA: 0 Median : 609.77 Median : -90.682 Mean : 714.79 Mean : -25.087 3rd Qu.: 912.27 3rd Qu.: 140.087 Max. :1969.57 Max. :1485.503 $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.15 Min. : 2251.5 POINT :3103 Min. : 47.450 1st Qu.: 111.52 1st Qu.: 4679.9 epsg:NA: 0 1st Qu.: 68.410 Median : 283.45 Median : 6688.4 Median : 81.783 Mean : 348.14 Mean : 8516.3 Mean : 88.646 3rd Qu.: 535.64 3rd Qu.:10854.7 3rd Qu.:104.186 Max. :1338.23 Max. :29401.2 Max. :171.468 exp_var: np dist gamma dir.hor dir.ver id 1 17 59.334696 23959.646 0 0 var1 2 36 86.014494 16671.647 0 0 var1 3 114 131.028697 33270.313 0 0 var1 4 149 176.188454 52846.333 0 0 var1 5 184 226.756517 34555.290 0 0 var1 6 711 337.603595 50150.725 0 0 var1 7 830 502.047734 50347.613 0 0 var1 8 1349 713.214854 61255.291 0 0 var1 9 1314 961.271787 70311.908 0 0 var1 10 1139 1213.411568 75870.335 0 0 var1 11 1355 1506.550518 72043.880 0 0 var1 var_model: model psill range kappa 1 Nug 2454.9391 0.00000 0.0 2 Ste 68746.1473 482.02453 0.3 Sums of squares betw. var. model and sample var.[1] 2044181.5 > pos = posPredictionInterval(kriging_result, 95, 75) > summary(pos) Upper and lower boundary and the position relative to value: upper lower position Min. : 23.995 Min. :-490.572 higher :1790 1st Qu.: 276.975 1st Qu.: -47.453 lower : 107 Median : 433.154 Median : 131.751 not distinguishable:1206 Mean : 521.886 Mean : 174.399 3rd Qu.: 732.911 3rd Qu.: 359.361 Max. :1480.025 Max. :1196.429 geometry POINT :3103 epsg:NA: 0 $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 22.59 0.31 22.89