[R-sig-Geo] optimizeNetwork function in intamapInteractive package (how to use auxiliary predictors)

[Tomislav Hengl]

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/