## Commands for descriptive spatial analysis using geoR
##
##   Updated for R Version 2.0 by Sudipto Banerjee, March 2008
##  

rm(list=ls())

library(akima)
library(MBA)
library(geoR)

##This is for pretty image maps, but no good for B&W
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", 
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

##Commands on p61

input.data <- read.table("http://www.biostat.umn.edu/~sudiptob/DukeSummerCourse/FORMGMT.txt", header=T)

input.data <- input.data[1:150,]

input.data$Easting <- (input.data$Longi*pi/180)*6371
input.data$Northing <- (input.data$Lat*pi/180)*6371

coords <- cbind(input.data$Easting, input.data$Northing)

interp.Y <- interp(input.data$Easting, input.data$Northing, input.data$Y)
image(interp.Y, xlim=range(input.data$Easting),ylim=range(input.data$Northing), col=jet.colors(200))

readline("\nHit enter to add contours to this image plot:\n")
contour(interp.Y, add=T)

mba.surf.Y <- mba.surf(cbind(coords, input.data$Y), no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(interp.Y, xlim=range(input.data$Easting),ylim=range(input.data$Northing), col=jet.colors(200))

readline("\nHit enter to add contours to this image plot:\n")
contour(interp.Y, add=T)

##readline("\nHit enter for perspective (3-d) plot:\n")
##persp(int.scp, xlim=range(myscallops$long),ylim=range(myscallops$lat),theta=0)

## Commands on p62

readline("\n\n\nHit enter for empirical semivariogram plotting:\n")

obj <- cbind(coords, input.data$Y)
Y.geo <- as.geodata(obj, coords.col=1:2, data.col=3)

Y.var <- variog(Y.geo, estimator.type="classical")

Y.var.robust <- variog(Y.geo, estimator.type="modulus")

par(mfrow=c(2,1))
plot(Y.var)
plot(Y.var.robust)

## Commands on p63

readline("\n\n\nHit enter for Matern variogram fitting:\n")
Y.var.fit <- variofit(Y.var.robust, ini.cov.pars=c(1.0,2.0),cov.model="matern",fix.nugget=FALSE, kappa=0.5, nugget=1.0)

print(summary(Y.var.fit))


readline("\n\n\nHit enter for Likelihood fitting:\n")
Y.lik.fit <- likfit(Y.geo, ini.cov.pars=c(1.0,2.0),cov.model="exponential",trend="cte", fix.nugget=FALSE,nugget=1.0, nospatial=FALSE)
print(summary(Y.lik.fit))


answer <- readline("Enter `1' for Bayesian kriging (TIME CONSUMING!):  ")
if (answer==1) {
scallops.bayes1 <- krige.bayes(Y.geo, locations="no", borders=NULL, model=model.control(trend.d = "cte", cov.model="matern", kappa=0.5, lambda=1), prior=prior.control(beta.prior="flat",sigmasq.prior="reciprocal",tausq.rel.prior="uniform",tausq.rel.discrete=seq(from=0.0,to=1.0,by=0.1), tausq.rel=0.5), output=output.control(n.posterior=2000, quantile=TRUE))

out <- scallops.bayes1$posterior ## This is wrong in the text! Change scallops.krige.bayes to scallops.bayes1 
out <- out$sample
beta.qnt <- quantile(out$beta, c(0.50,0.025,0.975))
phi.qnt <- quantile(out$phi, c(0.50,0.025,0.975))
sigmasq.qnt <- quantile(out$sigmasq, c(0.50,0.025,0.975))
tausq.rel.qnt <- quantile(out$tausq.rel, c(0.50,0.025,0.975))

cat("\nSummaries for beta:\n");  print(beta.qnt)
cat("\nSummaries for phi:\n");   print(phi.qnt)
cat("\nSummaries for sigmasq:\n");  print(sigmasq.qnt)
cat("\nSummaries for tausq.rel:\n");  print(tausq.rel.qnt)

}
cat("\n End of booknew.R computer tutorial. \n\n")