Gaussian cosimulation doesn't produce the correct correlated simulations

[edzer]

Example script (see Gerard Heuvelink's post on r-sig-geo):

library(gstat)
library(sp)
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)]))