if (FALSE) {
# Comparison of Co in March and September
   
with(mosses, {
	
   nbins <- 12
   vgm.m <- sm.variogram(loc.m, Co.m, nbins = nbins, original.scale = TRUE,
                        ylim = c(0, 1.5))
   vgm.s <- sm.variogram(loc.s, Co.s, nbins = nbins, original.scale = TRUE,
                        add = TRUE, col.points = "blue")
                        
   trns <- function(x) (x / 0.977741)^4
   del <- 1000
   plot(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean), type = "b",
         ylim = c(0, 1.5), xlab = "Distance", ylab = "Semi-variogram")
   points(vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean), type = "b",
         col = "blue", pch = 2, lty = 2)
   plot(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean), type = "b",
         ylim = c(0, 1.5), xlab = "Distance", ylab = "Semi-variogram")
   points(vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean), type = "b",
         col = "blue", pch = 2, lty = 2)
   segments(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean - 2 * vgm.m$se),
         vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean + 2 * vgm.m$se))
   segments(vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean - 2 * vgm.s$se),
         vgm.s$distance.mean - del, trns(vgm.s$sqrtdiff.mean + 2 * vgm.s$se),
         col = "blue", lty = 2)
   mn <- (vgm.m$sqrtdiff.mean + vgm.s$sqrtdiff.mean) / 2
   se <- sqrt(vgm.m$se^2 + vgm.s$se^2)
   plot(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean), type = "n",
        ylim = c(0, 1.5), xlab = "Distance", ylab = "Semi-variogram")
   polygon(c(vgm.m$distance.mean, rev(vgm.m$distance.mean)),
        c(trns(mn - se), rev(trns(mn + se))),
        border = NA, col = "lightblue")  
   points(vgm.m$distance.mean, trns(vgm.m$sqrtdiff.mean))
   points(vgm.s$distance.mean, trns(vgm.s$sqrtdiff.mean), col = "blue", pch = 2)
   vgm1 <- sm.variogram(loc.m, Co.m, nbins = nbins, varmat = TRUE, 
                        display = "none")
   vgm2 <- sm.variogram(loc.s, Co.s, nbins = nbins, varmat = TRUE,
                        display = "none")
   nbin  <- length(vgm1$distance.mean)
   vdiff <- vgm1$sqrtdiff.mean - vgm2$sqrtdiff.mean
   tstat <- c(vdiff %*% solve(vgm1$V + vgm2$V) %*% vdiff)
   pval  <- 1 - pchisq(tstat, nbin)
   print(pval)
})
# Assessing isotropy for Hg in March
with(mosses, {
   sm.variogram(loc.m, Hg.m, model = "isotropic")
})
# Assessing stationarity for Hg in September
with(mosses, {
   vgm.sty <- sm.variogram(loc.s, Hg.s, model = "stationary")
   i <- 1
   image(vgm.sty$eval.points[[1]], vgm.sty$eval.points[[2]], vgm.sty$estimate[ , , i],
         col = topo.colors(20))
   contour(vgm.sty$eval.points[[1]], vgm.sty$eval.points[[2]], vgm.sty$sdiff[ , , i],
         col = "red", add = TRUE)
})
}
Run the code above in your browser using DataLab