################################################### ### chunk number 1: libs ################################################### options(width = 67) library(spBayes) library(MBA) library(geoR) library(fields) library(cluster) 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$ALLBIO02_KGH>0,] bio <- BEF.dat$ALLBIO02_KGH*0.001; log.bio <- as.matrix(log(0.001*BEF.dat[,c("BOLE02_KGH", "BRANCH02_KGH", "FOLIAGE02_KGH")])) colnames(log.bio) <- c("log.bole.mt", "log.branch.mt", "log.foliage.mt") coords <- as.matrix(BEF.dat[,c("XUTM","YUTM")]) plot(coords, pch=19, cex=0.5, xlab="Easting (m)", ylab="Northing (m)") cov(log.bio) ################################################### ### chunk number 3: biomassLM ################################################### lm.bole <- lm(log.bio[,"log.bole.mt"] ~ ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat) lm.branch <- lm(log.bio[,"log.branch.mt"] ~ ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat) lm.foliage <- lm(log.bio[,"log.foliage.mt"] ~ ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, data=BEF.dat) summary(lm.bole) summary(lm.branch) summary(lm.foliage) ################################################### ### chunk number 4: knots ################################################### m <- 100 km.knots <- kmeans(coords, m)$centers cl.knots <- clara(coords, m)$medoids cd.knots <- cover.design(coords, nd=m)$design plot(coords, pch=19, cex=0.5, xlab="Easting (m)", ylab="Northing (m)") points(km.knots, pch=5, cex=1, col="blue") points(cl.knots, pch=6, cex=1, col="green") points(cd.knots, pch=7, cex=1, col="red") legend("bottomright", cex=1, pch=c(19,5,6,7), bty="n", col=c("black","blue","green","red"), legend=c("observations","kmeans","clara","cover.design")) ################################################### ### chunk number 5: spMvLM ################################################### q <- 3 ##Specify starting values and collect samples A.starting <- diag(0.2,q)[lower.tri(diag(1,q), TRUE)] L.starting <- diag(0.2,q)[lower.tri(diag(1,q), TRUE)] A.tuning <- rep(0.01,length(L.starting)) L.tuning <- rep(0.01,length(L.starting)) n.samples <- 10 model <- list(log.bio[,"log.bole.mt"]~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, log.bio[,"log.branch.mt"]~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3, log.bio[,"log.foliage.mt"]~ELEV+SLOPE+SUM_02_TC1+SUM_02_TC2+SUM_02_TC3 ) bef.spMvLM <- spMvLM(model, coords=coords, knots=km.knots, data=BEF.dat, starting=list("phi"=rep(3/500,q),"A"=A.starting,"L"=L.starting), sp.tuning=list("phi"=rep(0.3,q),"A"=A.tuning, "L"=L.tuning), priors=list("phi.Unif"=rep(c(3/2500,3/100),q), modified.pp=TRUE, "K.IW"=list(q+1, diag(0.1,q)), "Psi.IW"=list(q+1, diag(0.1,q))), cov.model="exponential", n.samples=n.samples, sub.samples=c(1,n.samples,1), verbose=TRUE, n.report=100) ##save(bef.spMvLM, file="R-data/bef.spMvLM") load(file="R-data/bef.spMvLM") p.samples <- bef.spMvLM$p.samples p.samples[,paste("phi_",1:q,sep="")] <- 3/p.samples[,paste("phi_",1:q,sep="")] colnames(p.samples)[colnames(p.samples) == paste("phi_",1:q,sep="")] <- paste("Eff. range ",1:q,sep="") round(summary(mcmc(p.samples))$quantiles[,c(1,3,5)],3) ################################################### ### chunk number 6: spLMWFitted ################################################### X <- bef.spMvLM$X beta <- t(p.samples[,1:bef.spMvLM$p]) w <- bef.spMvLM$sp.effects Y <- rowMeans(X%*%beta+w) bole <- Y[seq(1,length(Y),q)] branch <- Y[seq(2,length(Y),q)] foliage <- Y[seq(3,length(Y),q)] res <- 100 par(mfrow=c(2,2)) surf <- mba.surf(cbind(coords,bole), no.X=res, no.Y=res, extend=FALSE)$xyz.est image.plot(surf, main="Bole fitted values") points(coords) surf <- mba.surf(cbind(coords,branch), no.X=res, no.Y=res, extend=FALSE)$xyz.est image.plot(surf, main="Branch fitted values") points(coords) surf <- mba.surf(cbind(coords,foliage), no.X=res, no.Y=res, extend=FALSE)$xyz.est image.plot(surf, main="Foliage fitted values") points(coords) ################################################### ### chunk number 7: spPredict ################################################### x.range <- range(coords[,1]) y.range <- range(coords[,2]) 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"]], BEF.grids[["band3"]],BEF.grids[["band4"]],BEF.grids[["band5"]]) 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),] pred.X <- mkMvX(list(pred.covars, pred.covars, pred.covars)) bef.bio.pred <- spPredict(bef.spMvLM, start=1, end=2, pred.coords=pred.coords, pred.covars=pred.X, verbose=TRUE) ##save(bef.bio.pred, file="R-data/bef.bio.pred") load(file="R-data/bef.bio.pred") ################################################### ### chunk number 8: spPredMean ################################################### y.pred.mu <- apply(bef.bio.pred$y.pred, 1, mean) y.pred.sd <- apply(bef.bio.pred$y.pred, 1, sd) bole.pred.mu <- y.pred.mu[seq(1,length(y.pred.mu),q)] branch.pred.mu <- y.pred.mu[seq(2,length(y.pred.mu),q)] foliage.pred.mu <- y.pred.mu[seq(3,length(y.pred.mu),q)] bole.pred.sd <- y.pred.sd[seq(1,length(y.pred.sd),q)] branch.pred.sd <- y.pred.sd[seq(2,length(y.pred.sd),q)] foliage.pred.sd <- y.pred.sd[seq(3,length(y.pred.sd),q)] bio.pred.grid <- as.data.frame(list(x=pred.coords[,1], y=pred.coords[,2], bole.mu=bole.pred.mu, branch.mu=branch.pred.mu, foliage.mu=foliage.pred.mu, bole.sd=bole.pred.sd, branch.sd=branch.pred.sd, foliage.sd=foliage.pred.sd)) coordinates(bio.pred.grid) <- c("x", "y") # promote to SpatialPointsDataFrame gridded(bio.pred.grid) <- TRUE # promote to SpatialGridDataFrame toImage <- function(x){as.image.SpatialGridDataFrame(x)} res <- 100 par(mfrow=c(3,2)) surf <- mba.surf(cbind(coords,log.bio[,"log.bole.mt"]), no.X=res, no.Y=res, extend=FALSE)$xyz.est z.lim <- range(surf[["z"]], na.rm=TRUE) image.plot(surf, xaxs = "r", yaxs = "r", main="Observed bole biomass") image.plot(toImage(bio.pred.grid["bole.mu"]), zlim=z.lim, xaxs = "r", yaxs = "r", main="Predicted bole biomass") surf <- mba.surf(cbind(coords,log.bio[,"log.branch.mt"]), no.X=res, no.Y=res, extend=FALSE)$xyz.est z.lim <- range(surf[["z"]], na.rm=TRUE) image.plot(surf, xaxs = "r", yaxs = "r", main="Observed branch biomass") image.plot(toImage(bio.pred.grid["branch.mu"]), zlim=z.lim, xaxs = "r", yaxs = "r", main="Predicted branch biomass") surf <- mba.surf(cbind(coords,log.bio[,"log.foliage.mt"]), no.X=res, no.Y=res, extend=FALSE)$xyz.est z.lim <- range(surf[["z"]], na.rm=TRUE) image.plot(surf, xaxs = "r", yaxs = "r", main="Observed foliage biomass") image.plot(toImage(bio.pred.grid["foliage.mu"]), zlim=z.lim, xaxs = "r", yaxs = "r", main="Predicted foliage biomass")