[R-sig-Geo] [gstat] simple reboxing test fails

[Tom Roche]

summary: I'm trying to test reboxing, or 3D interpolation, using
package=gstat. The idea is to create a uniform input grid of box/voxel
dimensions=2x2x2 with each box value=1, then "rebox" to a uniform output
grid with dimensions=1x1x8; presumably the output box values should
approximate the input box values. (Or am I missing something?) But all
my predicted box values=NA; how to fix?

details:

sessionInfo()
# gets
#> R version 2.15.1 (2012-06-22)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=C                 LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base
#> 
#> loaded via a namespace (and not attached):
#> [1] tools_2.15.1

library(gstat)

# horizontal grid sequences for 4x4 grid with cell area=4
hor.from.in <- 0
hor.to.in   <- 8
hor.by.in   <- 2
# grid corners
corners.hor.in <- seq(from=hor.from.in, to=hor.to.in, by=hor.by.in)
#corners.hor.in
# grid midpoints
midpts.hor.in <-
  seq(from=(hor.from.in+(hor.by.in/2)),
      to=(hor.to.in-(hor.by.in/2)),
      by=hor.by.in)
midpts.hor.in
#> [1] 1 3 5 7

# vertical grid sequence for 4 layers with height=2
ver.from.in <- 0
ver.to.in   <- 8
ver.by.in   <- 2
corners.ver.in <- seq(from=ver.from.in, to=ver.to.in, by=ver.by.in)
#corners.ver.in
midpts.ver.in <-
  seq(from=(ver.from.in+(ver.by.in/2)),
      to=(ver.to.in-(ver.by.in/2)),
      by=ver.by.in)
midpts.ver.in
#> [1] 1 3 5 7

grid3D.in <-
  expand.grid(x=midpts.hor.in, y=midpts.hor.in, z=midpts.ver.in)
# more on grid3D.in (the input grid) below

val.in <- 1.0 # data value for each input box of volume=8
data3D.in <-
  data.frame(v=rep(val.in, times=length(row.names(grid3D.in))))
coordinates(data3D.in) <- grid3D.in
summary(data3D.in)
#> Object of class SpatialPointsDataFrame
#> Coordinates:
#>   min max
#> x   1   7
#> y   1   7
#> z   1   7
#> Is projected: NA 
#> proj4string : [NA]
#> Number of points: 64
#> Data attributes:
#>   Min. 1st Qu.  Median  Mean 3rd Qu.  Max. 
#>      1       1       1     1       1     1 

# so each datum=1 in the input data (data3D.in)

gridded(grid3D.in)=~x+y+z
summary(grid3D.in)
#> Object of class SpatialPixels
#> Coordinates:
#>   min max
#> x   0   8
#> y   0   8
#> z   0   8
#> Is projected: NA 
#> proj4string : [NA]
#> Number of points: 64
#> Grid attributes:
#>   cellcentre.offset cellsize cells.dim
#> x                 1        2         4
#> y                 1        2         4
#> z                 1        2         4

# so there are 64 box/voxels each of volume=8 in the input grid

# horizontal grid sequences for 8x8 grid with cell area=1
hor.from.out <- 0
hor.to.out   <- 8
hor.by.out   <- 1
corners.hor.out <- seq(from=hor.from.out, to=hor.to.out, by=hor.by.out)
corners.hor.out
midpts.hor.out <- seq(from=(hor.from.out+(hor.by.out/2)), to=(hor.to.out-(hor.by.out/2)), by=hor.by.out)
midpts.hor.out
#> [1] 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5

# vertical grid sequence for 1 layer with height=8
ver.from.out <- 0
ver.to.out   <- 8
ver.by.out   <- 8
corners.ver.out <- seq(from=ver.from.out, to=ver.to.out, by=ver.by.out)
corners.ver.out
midpts.ver.out <- seq(from=(ver.from.out+(ver.by.out/2)), to=(ver.to.out-(ver.by.out/2)), by=ver.by.out)
midpts.ver.out
# 4

# This fails to create
# 64 boxes each of volume=8 (but 1x1x8 instead of 2x2x2)
grid3D.out <- expand.grid(x=midpts.hor.out, y=midpts.hor.out, z=midpts.ver.out)
gridded(grid3D.out)=~x+y+z
summary(grid3D.out)
# gets
#> Object of class SpatialPixels
#> Coordinates:
#>   min max
#> x 0.0 8.0
#> y 0.0 8.0
#> z 3.5 4.5 ############################ NO! I want cellsize=8
#> Is projected: NA 
#> proj4string : [NA]
#> Number of points: 64
#> Grid attributes:
#>   cellcentre.offset cellsize cells.dim
#> x               0.5        1         8
#> y               0.5        1         8
#> z               4.0        1         1

# This appears to create
# 64 boxes each of volume=8 (but 1x1x8 instead of 2x2x2)
grid3D.out <-
  SpatialGrid(
    GridTopology(
      c(0.5,0.5,4.0), # centers
      c(1.0,1.0,8.0), # sizes
      c(8,  8,  1)),  # dim/number
    proj4string = CRS(as.character(NA)))
summary(grid3D.out)
# gets
#> Object of class SpatialGrid
#> Coordinates:
#>      min max
#> [1,]   0   8
#> [2,]   0   8
#> [3,]   0   8
#> Is projected: NA 
#> proj4string : [NA]
#> Grid attributes:
#>   cellcentre.offset cellsize cells.dim
#> 1               0.5        1         8
#> 2               0.5        1         8
#> 3               4.0        8         1

# so each box/voxel volume=8 in the output grid (grid3D.out), no?

# Here I fail to get the values for the output grid :-(
data3D.out <- krige(formula=v ~ 1, data3D.in, grid3D.out)
summary(data3D.out)
# gets
#> Object of class SpatialGridDataFrame
#> Coordinates:
#>      min max
#> [1,]   0   8
#> [2,]   0   8
#> [3,]   0   8
#> Is projected: NA 
#> proj4string : [NA]
#> Grid attributes:
#>   cellcentre.offset cellsize cells.dim
#> 1               0.5        1         8
#> 2               0.5        1         8
#> 3               4.0        8         1
#> Data attributes:
#>    var1.pred    var1.var  
#>  Min.   :1   Min.   : NA  
#>  1st Qu.:1   1st Qu.: NA  
#>  Median :1   Median : NA  
#>  Mean   :1   Mean   :NaN  
#>  3rd Qu.:1   3rd Qu.: NA  
#>  Max.   :1   Max.   : NA  
#>              NA's   :64   

What should I call to get the values for the output grid boxes?

Your assistance is appreciated, Tom Roche <Tom_Roche at pobox.com>