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)
> library(gstat)
> set.seed(13331)
> 
> x = runif(10)
> y = runif(10)
> z = rnorm(10)
> d = data.frame(x=x,y=y,z=z)
> bl = c(1,1)
> 
> nd = data.frame(x=.5, y=.5)
> 
> # single block:
> library(sp)
> coordinates(d) = ~x+y
> nd = SpatialPoints(nd)
> xx = krige(z~1, d, nd, model=vgm(1, "Exp", 1), block = bl,
+ 	set = list(nb=4))
[using ordinary kriging]
> 
> ring = cbind(c(0,bl[1],bl[1],0,0), c(0,0,bl[2],bl[2],0))
> r1 = SpatialPolygons(list(Polygons(list(Polygon(ring)), ID = "xx")))
> a = data.frame(a = 1, b = 2)
> rownames(a) = "xx"
> r1df = SpatialPolygonsDataFrame(r1, a)
> 
> g = gstat(formula=z~1, data=d, model=vgm(1, "Exp", 1))
> args = list(type = "regular", n=16, offset=c(0.5,0.5))
> yy = predict(g, r1df, block = bl, sps.args = args)
[using ordinary kriging]
> 
> xx = as.data.frame(xx)
> yy = as.data.frame(yy)
> row.names(xx) = row.names(yy)
> print(all.equal(xx,yy))
[1] TRUE
> 
> ## multiple blocks of equal size:
> nd = data.frame(x= 0:4 + .5, y=rep(.5,5))
> nd = SpatialPoints(nd)
> xx = krige(z~1, d, nd, model=vgm(1, "Exp", 1), block = bl,
+ 	set = list(nb=4))
[using ordinary kriging]
> ring0 = cbind(c(0,bl[1],bl[1],0,0), c(0,0,bl[2],bl[2],0))
> ring1 = cbind(c(1+0,1+bl[1],1+bl[1],1+0,1+0), c(0,0,bl[2],bl[2],0))
> ring2 = cbind(c(2+0,2+bl[1],2+bl[1],2+0,2+0), c(0,0,bl[2],bl[2],0))
> ring3 = cbind(c(3+0,3+bl[1],3+bl[1],3+0,3+0), c(0,0,bl[2],bl[2],0))
> ring4 = cbind(c(4+0,4+bl[1],4+bl[1],4+0,4+0), c(0,0,bl[2],bl[2],0))
> r1 = SpatialPolygons(list(
+ 	Polygons(list(Polygon(ring0)), ID = "x0"),
+ 	Polygons(list(Polygon(ring1)), ID = "x1"),
+ 	Polygons(list(Polygon(ring2)), ID = "x2"),
+ 	Polygons(list(Polygon(ring3)), ID = "x3"),
+ 	Polygons(list(Polygon(ring4)), ID = "x4")
+ 	))
> df = data.frame(a=rep(1,5), b= rep(1,5))
> rownames(df) = c("x0", "x1", "x2", "x3", "x4")
> r1df = SpatialPolygonsDataFrame(r1, df)
> 
> yy = predict(g, r1, block = bl, sps.args = args)
[using ordinary kriging]
> xx = as.data.frame(xx)
> yy = as.data.frame(yy)
> row.names(xx) = row.names(yy)
> all.equal(xx, yy)
[1] TRUE
> 
> ## multiple blocks of equal size:
> args = list(type = "regular", cellsize=.25, offset=c(0.5,0.5), n=16)
> yy = predict(g, r1, block = bl, sps.args = args)
[using ordinary kriging]
> xx = as.data.frame(xx)
> yy = as.data.frame(yy)
> row.names(xx) = row.names(yy)
> print(all.equal(xx, yy))
[1] TRUE
> 
> ## multiple blocks of varying size:
> nd = data.frame(x=c(0.5, 2, 4.5), y=c(0.5, 1, 1.5))
> nd = SpatialPoints(nd)
> bl = c(1,1)
> ring0 = cbind(c(0,bl[1],bl[1],0,0), c(0,0,bl[2],bl[2],0))
> xx1 = krige(z~1, d, nd[1], model=vgm(1, "Exp", 1), block = bl,
+ 	set = list(nb=4))
[using ordinary kriging]
> bl = c(2,2)
> ring1 = cbind(c(1+0,1+bl[1],1+bl[1],1+0,1+0), c(0,0,bl[2],bl[2],0))
> xx2 = krige(z~1, d, nd[2], model=vgm(1, "Exp", 1), block = bl,
+ 	set = list(nb=4))
[using ordinary kriging]
> bl = c(3,3)
> ring2 = cbind(c(3+0,3+bl[1],3+bl[1],3+0,3+0), c(0,0,bl[2],bl[2],0))
> xx3 = krige(z~1, d, nd[3], model=vgm(1, "Exp", 1), block = bl,
+ 	set = list(nb=4))
[using ordinary kriging]
> r1 = SpatialPolygons(list(
+ 	Polygons(list(Polygon(ring0)), ID = "x0"),
+ 	Polygons(list(Polygon(ring1)), ID = "x1"),
+ 	Polygons(list(Polygon(ring2)), ID = "x2")
+ 	))
> df = data.frame(a = rep(1,3), b = rep(1,3))
> rownames(df) = c("x0", "x1", "x2")
> r1df = SpatialPolygonsDataFrame(r1, df)
> 
> args = list(type = "regular", n=16, offset=c(0.5,0.5))
> yy = predict(g, r1df, block = bl, sps.args = args)
[using ordinary kriging]
> xx = rbind(as.data.frame(xx1), as.data.frame(xx2), as.data.frame(xx3))
> row.names(xx) = 1:3
> xx = as.data.frame(xx)
> yy = as.data.frame(yy)
> row.names(xx) = row.names(yy)
> print(all.equal(xx, yy))
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   1.09    0.28    1.34