## Commands in the geoR tutorial, Banerjee, Carlin, and Gelfand text, pp.61-65
##
##   Updated for R Version 2.0 by Sudipto Banerjee, March 2005
##   Monkeyed With by Brad Carlin, June 2005

set.seed(1234)
library(akima)
par(mfrow=c(1,1))

##Commands on p61

myscallops <- read.table("http://www.biostat.umn.edu/~brad/data/myscallops.txt", header=T)

int.scp <- interp(myscallops$long, myscallops$lat, myscallops$lgcatch)
image(int.scp, xlim=range(myscallops$long),ylim=range(myscallops$lat))

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

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

## Do both of these plots using both "interp" (linear interpolation) and
## "interp.new" (spline interpolation)

readline("\n\n\nHit enter for a comparison of these plots under 'interp' and 'interp.new':\n")

min.lgcatch<-min(myscallops$lgcatch)
dif.lgcatch<-max(myscallops$lgcatch)-min(myscallops$lgcatch)
int.scp.new <- interp.new(myscallops$long, myscallops$lat, myscallops$lgcatch)

par(mfrow=c(2,2))
zlim<-range(int.scp$z,int.scp.new$z,na.rm=T)
image(int.scp, xlim=range(myscallops$long),ylim=range(myscallops$lat),zlim=zlim,main="interp")
points(myscallops$long, myscallops$lat,pch=24,bg=hsv(0.6+(myscallops$lgcatch-min.lgcatch)/dif.lgcatch*0.4))
persp(int.scp, xlim=range(myscallops$long),ylim=range(myscallops$lat),theta=0,zlim=zlim,main="interp")

image(int.scp.new, xlim=range(myscallops$long),ylim=range(myscallops$lat),zlim=zlim,main="interp.new")
points(myscallops$long, myscallops$lat,pch=24,bg=hsv(0.6+(myscallops$lgcatch-min.lgcatch)/dif.lgcatch*0.4))
persp(int.scp.new, xlim=range(myscallops$long),ylim=range(myscallops$lat),theta=0,zlim=zlim,main="interp.new")


## Commands on p62

readline("\n\n\nHit enter for empirical semivariogram plotting:\n")
library(geoR) ##This seems to be missing.
R <- 6371
obj <- cbind(myscallops$long, myscallops$lat, myscallops$lgcatch)
scallops.geo <- as.geodata(obj, coords.col=1:2, data.col=3) ## This is wrong in the text! Change myscallops to obj 
scallops.var <- variog(scallops.geo, estimator.type="classical")

scallops.var ##This just prints out a variogram object. Am not sure whether we need this.

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

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

## Commands on p63

readline("\n\n\nHit enter for Matern variogram fitting:\n")
scallops.var.fit <- variofit(scallops.var.robust,ini.cov.pars=c(1.0,2.0),cov.model="matern",fix.nugget=FALSE,kappa=0.5,nugget=1.0) ## To be consistent with the output change the current version to this.

print(summary(scallops.var.fit))

##> scallops.var.fit
##variofit: model parameters estimated by WLS (weighted least squares):
##covariance model is: matern with fixed kappa = 0.5 (exponential)
##parameter estimates:
##  tausq sigmasq     phi 
## 0.0000  5.1110  0.2065 
## variofit: minimised weighted sum of squares = 4017.547


## Commands on p64

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

#   NOTE:  may need to avoid "summary" and simply type
#print(scallops.lik.fit)
#   with some versions of R


##> scallops.lik.fit 
## likfit: estimated model parameters:
##    beta    tausq  sigmasq      phi 
## "2.3748" "0.0947" "5.7675" "0.2338" 
## likfit: maximised log-likelihood = -285.9


## Commands on p65

answer <- readline("Enter `1' for Bayesian kriging (TIME CONSUMING!):  ")
if (answer==1) {
scallops.bayes1 <- krige.bayes(scallops.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.01)))

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