###################################################
### chunk number 1: libs
###################################################
options(width = 65) 
library(spBayes)
library(MBA)
library(geoR)
library(fields)
library(rgl)
library(sp)
library(maptools)
library(rgdal)
library(classInt)
library(gstat)
library(lattice)
library(xtable)


###################################################
### chunk number 2: coords
###################################################
data(BEF.dat)

BEF.dat <- BEF.dat[BEF.dat$BAREA02_TOT>0,]

EH <- as.integer(BEF.dat[,"EH_02BAREA"] > 0.25)

coords <- as.matrix(BEF.dat[,c("XUTM","YUTM")])

plot(coords, typ="n", xlab="Easting (m)", ylab="Northing (m)")
points(coords[as.logical(EH),], pch=19, cex=1)
points(coords[!as.logical(EH),], pch=1, cex=1)
legend("bottomright", pch=c(19, 1), cex=c(1,1), legend=c("High EH basal area", "Low EH basal area"), bty="n")


###################################################
### chunk number 3: EHlogit
###################################################
x.res <- 200; y.res <- 200

EH.logit <- glm(EH ~ ELEV+SLOPE, data=BEF.dat, family=binomial("logit"))
summary(EH.logit)



###################################################
### chunk number 4: spGLM
###################################################
beta.starting <- coefficients(EH.logit)
beta.tuning <- t(chol(vcov(EH.logit)))

n.samples <- 10
bef.starting <- spGLM(EH~ELEV+SLOPE, data=BEF.dat, coords=coords, 
                      starting=list("beta"=beta.starting, "phi"=3/600, "sigma.sq"=10, "w"=0),
                      tuning=list("beta"=beta.tuning, "phi"=0.0, "sigma.sq"=0.0, "w"=0.01),
                      priors=list("phi.Unif"=c(3/2000, 3/10), "sigma.sq.IG"=c(2, 10)),
                      cov.model="exponential",
                      n.samples=n.samples, verbose=TRUE, n.report=100)

burn.in <- as.integer(0.75*n.samples)
w.starting <- rowMeans(bef.starting$sp.effects[,burn.in:n.samples])

##save(w.starting, file = "R-data/w.starting")
load(file = "R-data/w.starting")

bef.sp <- spGLM(EH~ELEV+SLOPE, data=BEF.dat, coords=coords, 
                starting=list("beta"=beta.starting, "phi"=3/600, "sigma.sq"=5, "w"=w.starting),
                tuning=list("beta"=beta.tuning, "phi"=0.05, "sigma.sq"=0.03, "w"=0.01),
                priors=list("phi.Unif"=c(3/1000, 3/10), "sigma.sq.IG"=c(2, 5)),
                cov.model="exponential",
                n.samples=n.samples, verbose=TRUE, n.report=500)

##save(bef.sp, file = "R-data/bef.sp")
load(file = "R-data/bef.sp")

n.samples <- nrow(bef.sp$p.samples)

burn.in <- as.integer(0.50*n.samples)
summary(mcmc(bef.sp$p.samples))$quantiles


###################################################
### chunk number 5: EHFittedSurf
###################################################
logit.fitted <- function(eta){1/(1+exp(-eta))}

fitted <- apply(bef.sp$X%*%t(bef.sp$p.samples[burn.in:n.samples,1:3])+bef.sp$sp.effects[,burn.in:n.samples]
                ,1,logit.fitted)

fitted.mu <- apply(fitted,2,mean)

par(mfrow=c(1,2))
surf <- mba.surf(cbind(coords, EH), no.X=x.res, no.Y=y.res, extend=FALSE)$xyz.est
image.plot(surf, xaxs = "r", yaxs = "r", main="Observed EH")

surf <- mba.surf(cbind(coords, fitted.mu), no.X=x.res, no.Y=y.res, extend=FALSE)$xyz.est
image.plot(surf, xaxs = "r", yaxs = "r", main="Fitted probabilities")


###################################################
### chunk number 6: spPredict
###################################################
BEF.shp <- readShapePoly("BEF-data/BEF_bound.shp")
BEF.poly <- as.matrix(BEF.shp@polygons[[1]]@Polygons[[1]]@coords)

BEF.grids <- readGDAL("BEF-data/dem_slope_lolosptc_clip_60.img")

pred.covars <- cbind(BEF.grids[["band1"]],BEF.grids[["band2"]])

pred.covars <- cbind(rep(1, nrow(pred.covars)), pred.covars)
pred.coords <- SpatialPoints(BEF.grids)@coords

pred.covars <- pred.covars[pointsInPoly(BEF.poly, pred.coords),]
pred.coords <- pred.coords[pointsInPoly(BEF.poly, pred.coords),]

burn.in <- as.integer(0.5*n.samples)
bef.EH.pred <- spPredict(bef.sp, pred.coords=pred.coords, pred.covars=pred.covars, start=burn.in, end=burn.in+1, verbose=FALSE)

##save(bef.EH.pred, file = "R-data/bef.EH.pred")
load(file = "R-data/bef.EH.pred")


###################################################
### chunk number 7: EHspPredMean
###################################################
EH.mu <- apply(bef.EH.pred$y.pred, 1, mean)
EH.sd <- apply(bef.EH.pred$y.pred, 1, sd)

EH.w.mu <- apply(bef.EH.pred$w.pred, 1, mean)
EH.w.sd <- apply(bef.EH.pred$w.pred, 1, sd)

EH.pred.grid <- as.data.frame(list(x=pred.coords[,1], y=pred.coords[,2], EH.mu=EH.mu, EH.sd=EH.sd, EH.w.mu=EH.w.mu, EH.w.sd=EH.w.sd))

coordinates(EH.pred.grid) <- c("x", "y") # promote to SpatialPointsDataFrame
gridded(EH.pred.grid) <- TRUE            # promote to SpatialGridDataFrame

toImage <- function(x){as.image.SpatialGridDataFrame(x)}

par(mfrow=c(1,2))

plot(coords, typ="n", xlab="Easting (m)", ylab="Northing (m)", main="Observed EH inventory plots")
points(coords[as.logical(EH),], pch=19, cex=1)
points(coords[!as.logical(EH),], pch=1, cex=1)
legend("bottomright", pch=c(19, 1), cex=c(1,1), legend=c("High EH basal area", "Low EH basal area"), bty="n")

image.plot(toImage(EH.pred.grid["EH.mu"]), xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)", main="Mean predicted probability of EH")
points(coords)


###################################################
### chunk number 8: EHspPredSD
###################################################

par(mfrow=c(1,2))
image.plot(toImage(EH.pred.grid["EH.sd"]), xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)", main="SD predicted probability of EH")

image.plot(toImage(EH.pred.grid["EH.w.mu"]), xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)", main="Mean predicted spatial effect of EH")