[R-sig-Geo] multivariate sequential gaussian cosimulation using gstat in R

[William Savran]

Dear All,

I am new to gstat, R, and fairly new to applying geostatistical methods, so I have a few questions.  I am looking to simulate a set of 3 (out of 4) correlated stochastic fields, which are all related through a linear model of coregionalization. One of the fields will be used as conditional data for the simulation.  With that being said, if anything seems alarming to you please let me know. Also, if you have any tips/trick that would be helpful for a novice please let me know.

Right now, my algorithm is working as follows (almost exactly following the demo):

# define gstat object
sim.g <- gstat(id='slip', formula=slip.sc~1, data=sim, nmax = 30, maxdist = 200, beta=0, set = list(nocheck = 1))
sim.g <- gstat(sim.g, 'psv', psv.sc~1, sim, nmax = 30, maxdist = 200, beta=0)
sim.g <- gstat(sim.g, 'vrup', vrup.sc~1, sim, nmax = 30, maxdist = 200, beta=0)
sim.g <- gstat(sim.g, 'mu0', mu0.sc~1, sim, nmax = 30, maxdist = 200, beta=0)
sim.g <- gstat(sim.g, model=vgm(1,"Exp",1000,1,anis=c(90,0.5)), fill.all=T)

# fit lmc
sim.fit = fit.lmc(var, sim.g, correct.diagonal = 1.01)

# perform SGSim
z <- predict(sim.fit, newdata=xy, nsim=1, debug.level = -1)  If I specify fit.ranges = TRUE and fit.LMC = TRUE, I get the warning message 

"Warning messages:
1: In fit.variogram(object, model, fit.sills = fit.sills, fit.ranges = fit.ranges,  :
 No convergence after 200 iterations: try different initial values?
2: In predict.gstat(sim.fit, newdata = xy, nsim = 1, debug.level = -1) :
 No Intrinsic Correlation or Linear Model of Coregionalization found
Reason: ranges differ”

and a simulation still happens.  What is going on behind the scenes here when there is no acceptable model of coregionalization or intrinsic correlation?  Do I only want to simulate using a constant range for all my 

Main Question:
----------------------------
In the SGSim step, how can I apply pre-existing data which can be used as conditional in the simulation process?  Based on the ids present in the gstat object,  I will have a data set of mu0, and conditional on this I want to simulate “psv”, “vrup”, and “slip”.  I have tried putting the data in the sim.fit object as well as in newdata and it doesn’t seem to make a difference.  What am I doing wrong?

Thanks in advance!

Best,
- William Savran