a <- read.table("http://www.stat.ucla.edu/~nchristo/statistics_c173_c273/soil.txt", header=TRUE)

image.orig <- image
library(gstat)

c1 <- ifelse(a$lead > 100, 1 , 0)


q <- gstat(id="c1", formula=c1~1, locations=~x+y, data=a)


plot(variogram(q))

v.fit <- fit.variogram(variogram(q), vgm(0.2,"Sph", 1000, 0.04))

plot(variogram(q),v.fit)


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)) 


pr_ok <- krige(id="c1",c1~1, locations=~x+y, model=v.fit, 
data=a, newdata=grd)

#Set values < 0 equal to zero, and values > 1 equal to 1:
for(i in 1:length(pr_ok$c1.pred) ){if(pr_ok$c1.pred[i] <0) {pr_ok$c1.pred[i]=0}}

for(i in 1:length(pr_ok$c1.pred) ){if(pr_ok$c1.pred[i] >1) {pr_ok$c1.pred[i]=1}}
 

#Collapse the predicted values into a matrix: 
qqq <- matrix(pr_ok$c1.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))) 

#Raster map:
image.orig(seq(from=x.range[1], to=x.range[2], by=50), 
seq(from=y.range[1],to=y.range[2], by=50), qqq, 
xlab="Westto East",ylab="South to North", main="Predicted values") 

#Add the observed data points on the raster map:
points(a)



#We want to assign "NA" values for the region outside the observed data #points:
X <- seq(from=x.range[1], to=x.range[2], by=50)
Y <- seq(from=y.range[1], to=y.range[2], by=50)

dz.mask <- function(grid, pts, x.dom, y.dom, window, mitre=2) {
	N <- length(pts[ ,1]) ; mask <- array( NA, dim(grid) )
	for(j in 1:N) {
		dx <- abs(x.dom - pts$x[j])
		dy <- abs(y.dom - pts$y[j])
		d.Mx <- tcrossprod( dx , rep(1, length(dy)) )^mitre + 
		tcrossprod( rep(1, length(dx)), dy )^mitre
		mask[ d.Mx < window^mitre ] <- FALSE
	}
	return(mask+grid)
}

qqq.masked <- dz.mask(qqq, a, X, Y, 250)

image.orig(X,Y,qqq.masked)
points(a, cex=a$lead/mean(a$lead), pch=19, col="green")


#Add contour lines:
contour(seq(from=x.range[1], to=x.range[2], by=50), 
+ seq(from=y.range[1],to=y.range[2], by=50), qqq.masked, 
levels=seq(0, 1, by=0.1), add=TRUE, col="black", labcex=1)