# Cokriging, zero distance semivariance [gstat]

Trying cokriging with simulated data, I faced a problem that did not seem one in the demo(cokriging) with the meuse dataset:

I can't use fit.lmc() default but need to specify fit.method=1 as the default fit.method=7 refuses to work with "zero distance semivariance". This is not a trouble in itself, but I surely do something wrong and/or do not understand what I am doing... although I try to follow the textbooks (gstat vignette, other good ref).

My question: my understanding of cokriging is that our interpolation of a first kind of measurement is helped by a denser, second kind of measurement. Provided that they are both present in some locations, so that regression can be done. Why then, are zero distance cross-semivariance a problem? Shouldn't it be a necessity?

Also, in textbooks, the meuse dataset is used without this error but measurements of the different covariates always overlap (all locations include measurements of all covariates). How then can there not be zero distance cross-semivariance?

I do hope the MWE below illustrates my misunderstanding and that someone can clarify this to me!

set.seed(1)
data <- data.frame(x=runif(100), y=runif(100))
data$z1=c(rnorm(20, mean=data$x*data$y), rep(NA, 80)) data$z2=rnorm(100, mean=data$x*data$y)
cor(data$z1[1:20], data$z2[1:20])
coordinates(data) <- ~ x + y

vgm1 <- variogram(z1~1, data[1:10,])
fit1 <- fit.variogram(vgm1, vgm(psill=2, model = "Sph",
range=.2, nugget=0))
plot(vgm1, fit1)
vgm2 <- variogram(z2~1, data)
fit2 <- fit.variogram(vgm2, vgm(psill=2, model = "Sph",
range=.2, nugget=0))
plot(vgm2, fit2)

(g <- gstat(NULL, id = "z1", form = z1 ~ 1, data=data[1:20,], nmax=5))
(g <- gstat(g, id = "z2", form = z2 ~ 1, data=data, nmax=10))
v.cross <- variogram(g)
plot(v.cross, pl=T)
g <- gstat(g, fill.all=T, model=vgm(psill=2, model = "Sph",
range=.2, nugget=0))
(g <- fit.lmc(v.cross, g, fit.method=1, correct.diagonal = 1.01)) # works
(g <- fit.lmc(v.cross, g, correct.diagonal = 1.01)) # does not work
plot(variogram(g), model=g\$model)