#Fit a squared exponential covariance function
 fit <- variofit(parana.bin,c(250,5000),"gaussian")
 lines.variomodel(fit)

#Set up prediction grid and krige
 r1 <- range(parana$coords[,1])
 r2 <- range(parana$coords[,2])
 loci <- expand.grid(seq(0,1,l=21),seq(0,1,l=21))
 loci <- cbind(r1[1]+(r1[2]-r1[1])*loci[,1],r2[1]+(r2[2]-r2[1])*loci[,2])
 parana.control <- krige.control(obj.model=fit)
 parana.krige <-krige.conv(parana,locations=loci,krige=parana.control)
 image(parana.krige,loc=loci,col=gray(seq(1,0.1,l=20)))