library(gstat)
data(meuse)
coordinates(meuse) <- ~x+y
data(meuse.grid)
coordinates(meuse.grid) <- ~x+y

# code from ASDAR book:
v <- variogram(log(zinc) ~ 1, meuse)
v.fit <- fit.variogram(v, vgm(1, "Sph", 800, 1))

g.cc <- gstat(NULL, "log.zinc", log(zinc)~1, meuse, model = v.fit)
meuse.grid$distn <- meuse.grid$dist - mean(meuse.grid$dist) +
  mean(log(meuse$zinc))
vd.fit <- v.fit
vov <- var(meuse.grid$distn) / var(log(meuse$zinc))
vd.fit$psill <- v.fit$psill * vov
g.cc <- gstat(g.cc, "distn", distn ~ 1, meuse.grid, nmax = 1,
              model=vd.fit, 
              merge = c("log.zinc","distn"))
vx.fit <- v.fit
vx.fit$psill <- sqrt(v.fit$psill * vd.fit$psill) *
  cor(meuse$dist, log(meuse$zinc)) #$
g.cc <- gstat(g.cc, c("log.zinc", "distn"), model = vx.fit) 
x <- predict(g.cc, meuse.grid)

# modified to generate error:
g.cc <- gstat(NULL, "log.zinc", log(zinc)~1, meuse, model = v.fit,
              maxdist=50) # => some points have no value of the
                          # primary variable in the neighbourhood:
                          # seems to cause crash
meuse.grid$distn <- meuse.grid$dist - mean(meuse.grid$dist) +
  mean(log(meuse$zinc))
vd.fit <- v.fit
vov <- var(meuse.grid$distn) / var(log(meuse$zinc))
vd.fit$psill <- v.fit$psill * vov
g.cc <- gstat(g.cc, "distn", distn ~ 1, meuse.grid, nmax = 1,
              model=vd.fit, 
              merge = c("log.zinc","distn"))
vx.fit <- v.fit
vx.fit$psill <- sqrt(v.fit$psill * vd.fit$psill) *
  cor(meuse$dist, log(meuse$zinc)) #$
g.cc <- gstat(g.cc, c("log.zinc", "distn"), model = vx.fit) 
x <- predict(g.cc, meuse.grid)

# is this the case for local ordinary kriging too?
# Code from ASDAR book:
lz.ok <- krige(log(zinc)~1, meuse, meuse.grid, v.fit)
# modified:
lz.ok <- krige(log(zinc)~1, meuse, meuse.grid, v.fit, maxdist=40)
# => works nonetheless

# ordinary cokriging?
# book
g <- gstat(NULL, "logCd", log(cadmium)~1, meuse)
g <- gstat(g, "logCu", log(copper)~1, meuse)
g <- gstat(g, "logPb", log(lead)~1, meuse)
g <- gstat(g, "logZn", log(zinc)~1, meuse)
vm <- variogram(g)
vm.fit <- fit.lmc(vm, g, vgm(1, "Sph", 800, 1))

cok.maps <- predict(vm.fit, meuse.grid)

# modified
g <- gstat(NULL, "logCd", log(cadmium)~1, meuse, maxdist=40)
g <- gstat(g, "logCu", log(copper)~1, meuse, maxdist=40)
g <- gstat(g, "logPb", log(lead)~1, meuse, maxdist=40)
g <- gstat(g, "logZn", log(zinc)~1, meuse, maxdist=40)
vm <- variogram(g)
vm.fit <- fit.lmc(vm, g, vgm(1, "Sph", 800, 1))

cok.maps <- predict(vm.fit, meuse.grid)
