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

[William Savran]

>On 06/06/16 19:03, William Savran wrote:
>> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

> Your script is quite messy and incomplete (where do sim, var, and xy
> come from? Does the "If I specify..." comment really relate to the
> predict call? Which command generates which warning?), and is not
> reproducible.
>
> The way you specify conditioning data is correct, so it is the basis for
> conditioning your simulations. The warnings should not be ignored,
> though: have you tried plotting the variograms & cross variograms with
> the (wrongly?) fitted model?

> When you give this little insight in what you do, I can only reply with
> wild guesses.
> --
> Edzer Pebesma
> Institute for Geoinformatics  (ifgi),  University of M?nster
> Heisenbergstra?e 2, 48149 M?nster, Germany; +49 251 83 33081
> Journal of Statistical Software:   http://www.jstatsoft.org/
> Computers & Geosciences:   http://elsevier.com/locate/cageo/
> Spatial Statistics Society http://www.spatialstatistics.info

Hi Edzer,

I apologize about not asking my question clearly, however it appears your
wild guess was spot on.  Now my initial simulations seem to be working as
expected, conditioning on mu0, and simulating slip,psv,and vrup from
scratch. Right now I am creating a new gstat object and copying over the
results from the output of fit.lmc as follows:

# create new gstat object for simulation from output of fit.lmc
sim.d <- gstat(id='mu0', formula=mu0.sc~1, model=sim.fit$model$mu0,
data=sim, nmax=30, beta=0.0)
sim.d <- gstat(sim.d, id='slip', formula=slip.dummy~1,
model=sim.fit$model$slip,  dummy=TRUE, nmax=30, beta=0.0)
sim.d <- gstat(sim.d, id='psv', formula=psv.dummy~1,
model=sim.fit$model$psv, dummy=TRUE, nmax=30, beta=0.0)
sim.d <- gstat(sim.d, id='vrup', formula=vrup.dummy~1,
model=sim.fit$model$vrup, dummy=TRUE, nmax=30, beta=0.0)

# enter cross-variograms
sim.d <- gstat(sim.d, id=c("slip","psv"), model=sim.fit$model$slip.psv)
sim.d <- gstat(sim.d, id=c("slip","vrup"), model=sim.fit$model$slip.vrup)
sim.d <- gstat(sim.d, id=c("slip","mu0"), model=sim.fit$model$slip.mu0)
sim.d <- gstat(sim.d, id=c("psv","vrup"), model=sim.fit$model$psv.vrup)
sim.d <- gstat(sim.d, id=c("psv","mu0"), model=sim.fit$model$psv.mu0)
sim.d <- gstat(sim.d, id=c("vrup","mu0"), model=sim.fit$model$vrup.mu0)

As far as my question regarding the warnings from fit.lmc and predict is
concerned, I increased the complexity of the model slightly and include
variogram basis functions with different correlation lengths.   When I do
this I can achieve a better fit to the empirical variograms, and it
provides an acceptable linear model of coregionalization.  Thanks again for
your help!

Best regards,
William

	[[alternative HTML version deleted]]