R Under development (unstable) (2024-09-04 r87094 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 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.

> library(sp)
> data(meuse)
> coordinates(meuse) = c("x", "y")
> new.locs <- SpatialPoints(data.frame(
+ 	x = c(181170, 180310, 180205, 178673, 178770, 178270),
+ 	y = c(333250, 332189, 331707, 330066, 330675, 331075)))
> library(gstat)
> krige(zinc ~ 1,  meuse, new.locs, vgm(1.34e5, "Sph", 800, nug = 2.42e4), 
+ 		block = c(40,40), nmax = 40)
[using ordinary kriging]
       coordinates var1.pred  var1.var
1 (181170, 333250)  268.7576  19447.67
2 (180310, 332189)  663.4915  16991.33
3 (180205, 331707)  251.4606  21579.71
4 (178673, 330066)  532.5776  73807.91
5 (178770, 330675)  664.4039  23589.17
6 (178270, 331075)  565.5436 155113.23
> 
> new.locs <- SpatialPoints(data.frame(x = c(181170), y = c(333250)))
> 
> disc = c(-15,-5,5,15)
> block.irreg <- data.frame(expand.grid(x = disc, y = disc))
> block.irreg
     x   y
1  -15 -15
2   -5 -15
3    5 -15
4   15 -15
5  -15  -5
6   -5  -5
7    5  -5
8   15  -5
9  -15   5
10  -5   5
11   5   5
12  15   5
13 -15  15
14  -5  15
15   5  15
16  15  15
> 
> # first disable default Gaussian quadrature used for block integration, by
> # setting nblockdiscr explicitly:
> krige(zinc ~ 1, meuse, new.locs, model = vgm(1.34e5, "Sph", 800, nug = 2.42e4), 
+ 		block = c(40,40), nmax = 40, set = list(nblockdiscr=4))
[using ordinary kriging]
       coordinates var1.pred var1.var
1 (181170, 333250)  268.7324 19568.76
> # now pass the same block discretization as block.irreg:
> krige(zinc ~ 1, meuse, new.locs, vgm(1.34e5, "Sph", 800, nug = 2.42e4), 
+ 		block = block.irreg, nmax = 40)
[using ordinary kriging]
       coordinates var1.pred var1.var
1 (181170, 333250)  268.7324 19568.76
> # check weights argument:
> block.irreg <- data.frame(expand.grid(x = disc, y = disc), weights = rep(1/16, 16))
> krige(zinc ~ 1, meuse, new.locs, vgm(1.34e5, "Sph", 800, nug = 2.42e4), 
+ 		block = block.irreg, nmax = 40)
[using ordinary kriging]
       coordinates var1.pred var1.var
1 (181170, 333250)  268.7324 19568.76
> 
> proc.time()
   user  system elapsed 
   1.10    0.12    1.15