#Cokriging with one co-located variable: a <- read.table("http://www.stat.ucla.edu/~nchristo/statistics_c173_c273/soil.txt", header=TRUE) library(gstat) #Begin with the target variable log_lead: #Create a gstat object using the target variable first: g1 <- gstat(id="log_lead", formula = log(lead)~1, locations = ~x+y, data = a) #Append to g1 the co-located variable: g1 <- gstat(g1,id="log_zinc", formula = log(zinc)~1, locations = ~x+y, data = a) #Plot auto- and cross- variograms: plot(variogram(g1)) #==================================================¯ #Use only target variable to model the variogram: g <- gstat(id="log_lead", formula = log(lead)~1, locations = ~x+y, data = a) #Plot the variogram of the target variable: plot(variogram(g)) #Fit a model to the sample variogram: v.fit <- fit.variogram(variogram(g), vgm(0.5,"Sph",1000,0.1)) #Plot the model variogram: plot(variogram(g),v.fit) #==================================================¯ #Now go back the gstat object with the two variables (target and co-located): vm <- variogram(g1) #Fit model variogram to all the sample variograms. Important: v.fit is the fitted model we used above! vm.fit <- fit.lmc(vm, g1, model=v.fit) #Plot the fitted variograms to all the sample variograms: plot(vm, vm.fit) #Cokriging predictions: x.range <- as.integer(range(a[,1])) y.range <- as.integer(range(a[,2])) grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=50), y=seq(from=y.range[1], to=y.range[2], by=50)) ck <- predict(vm.fit, grd) names(ck) ck$log_lead.pred #Plots: #Collapse the predicted values into a matrix: qqq <- matrix(ck$log_lead.pred, length(seq(from=x.range[1], to=x.range[2], by=50)), length(seq(from=y.range[1], to=y.range[2], by=50))) image(seq(from=x.range[1], to=x.range[2], by=50), seq(from=y.range[1], to=y.range[2], by=50), qqq) #Add contours: #Find first the range of the values: range(ck$log_lead.pred) contour(seq(from=x.range[1], to=x.range[2], by=50), seq(from=y.range[1],to=y.range[2], by=50), qqq, add=TRUE, col="black", levels=seq(3.5, 6.3, by=0.2), labcex=1) points(a)