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

[Jon Olav Skoien]

Hi Tom,

it is great that you like the intamap-packages!
The optimization function could not use auxiliary predictors properly in 
the version you tested, but there is now a new version on CRAN (it might 
take some time before it is distributed) where it is possible to do what 
you ask for.

Please note that there is an approximation in the function for 
estimating the optimization criteria. It is in this implementation 
necessary to know the auxiliary variables at the candidate locations, 
and these are estimated from the prediction grid. This should for most 
applications only have a minor effect.

There is an example in the documentation of optimizeNetwork.

Cheers,
Jon

Tomislav Hengl wrote:
>
> 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/
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo