optimizeNetwork function in intamapInteractive package (how to use auxiliary predictors)

Dear Edzer, Jon, Stephanie, Olivier et al.,

First, let me congratulate you on your new package intamap / 
intamapInteractive that I feel will significantly advance practice of 
geostatistical mapping. Really a great contribution (I knew about SSA 
already in 2000, but was always lacking an operational tool).

My interest at the moment are the sampling optimization algorithms and 
their use in combination with universal kriging / regression-kriging. 
I had no problems to run spatially simulated annealing with your 
package (OK model), but could not grasp how to do the same thing using 
a universal kriging model (auxiliary predictors). Here is the code:

library(gstat)
library(RSAGA)
library(maptools)
library(intamapInteractive)

# load data:

data(meuse)
coordinates(meuse) <- ~x+y
data(meuse.grid)
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE

# estimate variograms (OK/UK):

vt.fit <- fit.variogram(variogram(log1p(zinc)~1, meuse), vgm(1, 
"Exp", 300, 1))
vr.fit <- fit.variogram(variogram(log1p(zinc)~sqrt(dist), meuse), 
vgm(1, "Exp", 300, 1))

# study area of interest:

write.asciigrid(meuse.grid["mask"], "meuse_mask.asc")
rsaga.esri.to.sgrd(in.grids="meuse_mask.asc", 
out.sgrd="meuse_mask.sgrd", in.path=getwd())

# raster to polygon conversion;

rsaga.geoprocessor(lib="shapes_grid", module=6, 
param=list(GRID="meuse_mask.sgrd", SHAPES="meuse_mask.shp", 
CLASS_ALL=1))
mask <- readShapePoly("meuse_mask.shp", 
proj4string=CRS(as.character(NA)), force_ring=T)

# prepare observations / predGrid objects:

observations <- meuse["zinc"]
attr(observations at coords, "dimnames")[[2]] <- c("x", "y")
predGrid <- data.frame(meuse.grid["dist"])
predGrid$dist <- NULL  # predGrid has to be "SpatialPoints" type?
coordinates(predGrid) <- ~x+y

windows()
# add 20 more points assuming OK model (SSA method):

optim.ok <- optimizeNetwork(observations, predGrid, candidates=mask, 
method="ssa", action="add", nDiff=20, model=vt.fit, criterion="MUKV", 
plot=FALSE, nr_iterations=50, nmax=40)

This works fine. But how do I now optimize samples using universal 
kriging model (with auxiliary predictors e.g. "dist")? I simply could 
not figure out from your examples. I also tried looking at the Package 
source R functions, but could not get it running.

It should be something like this:

observations <- pts.ov[c("zinc", "dist")]
predGrid <- data.frame(meuse.grid["dist"])
coordinates(predGrid) <- ~x+y
optim.rk <- optimizeNetwork(observations, predGrid, candidates=mask, 
method="ssa", action="add", nDiff=20, model=vr.fit, criterion="MUKV", 
plot=FALSE, nr_iterations=50, formula=log1p(zinc)~sqrt(dist), nmax=40)

thanks!

Tom Hengl
http://home.medewerker.uva.nl/t.hengl/