[R-sig-Geo] cokriging simulation

[Heuvelink, Gerard]

Dear List,

We wish to create simulations from a cokriging model. I used to be able to do that in the past, but now I am struggling to get the correlations preserved. I dimplified the problem to a simple case with cross-correlation 0.8 (see script below), but simulations show no correlation at all. What do I do wrong? Thank you for your help, Gerard

library(gstat)
set.seed(12345)

d <- data.frame(0,0,0,0)
names(d) <- c("x","y","lmCres","lmNres")
coordinates(d) = ~x+y

nugC <- 1; nugN <- 1; nugCN <- 0.8
psillC <- 5; psillN <- 5; psillCN <- 4
range <- 10

g_l <- gstat(id=c("lmCres"),formula=lmCres~1,data=d,beta=0,nmax=100,
             model=vgm(psillC,"Sph",range=10,nugC))
g_l <- gstat(g_l,id="lmNres",formula=lmNres~1,data=d,beta=0,nmax=100,
             model=vgm(psillN,"Sph",range=10,nugN))
g_l <- gstat(g_l,id=c("lmCres","lmNres"),model=vgm(psillCN,"Sph",range=10,nugCN))

df <- data.frame(x=1,y=1)
coordinates(df) <- ~x+y

nsim <- 250
sim <- predict(g_l,df,nsim=nsim)
simdf <- as.data.frame(sim)

plot(as.numeric(simdf[1,3:(nsim+2)]),as.numeric(simdf[1,(nsim+3):((2*nsim)+2)]))
cor(as.numeric(simdf[1,3:(nsim+2)]),as.numeric(simdf[1,(nsim+3):((2*nsim)+2)]))

	[[alternative HTML version deleted]]