################################################### ### 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) ################################################### ### chunk number 2: coords ################################################### data(WEF.dat) WEF.dat <- WEF.dat[!apply(WEF.dat[,c("East_m","North_m","DBH_cm","Tree_height_m","ELEV_m")], 1, function(x)any(is.na(x))),] DBH <- WEF.dat$DBH_cm HT <- WEF.dat$Tree_height_m coords <- as.matrix(WEF.dat[,c("East_m","North_m")]) plot(coords, pch=1, cex=sqrt(DBH)/10, col="darkgreen", xlab="Easting (m)", ylab="Northing (m)") leg.vals <- round(quantile(DBH),0) legend("topleft", pch=1, legend=leg.vals, col="darkgreen", pt.cex=sqrt(leg.vals)/10, bty="n", title="DBH (cm)") ################################################### ### chunk number 3: colorRamp ################################################### col.br <- colorRampPalette(c("blue", "cyan", "yellow", "red")) col.pal <- col.br(5) ################################################### ### chunk number 4: classInt ################################################### quant <- classIntervals(DBH, n=5, style="quantile") fisher <- classIntervals(DBH, n=5, style="fisher") kmeans <- classIntervals(DBH, n=5, style="kmeans") bclust <- classIntervals(DBH, n=5, style="bclust") par(mfrow=c(2,2)) plot(quant, pal=col.pal, xlab="DBH", main="Quantile") plot(fisher, pal=col.pal, xlab="DBH", main="Fisher") plot(kmeans, pal=col.pal, xlab="DBH", main="Kmeans") plot(bclust, pal=col.pal, xlab="DBH", main="Bclust") ################################################### ### chunk number 5: quantFisher ################################################### quant.col <- findColours(quant, col.pal) fisher.col <- findColours(fisher, col.pal) par(mfrow=c(1,2)) plot(coords, col=quant.col, pch=19, cex=0.5, main="Quantile", xlab="Easting (m)", ylab="Northing (m)") legend("topleft", fill=attr(quant.col, "palette"), legend=names(attr(quant.col, "table")), bty="n") plot(coords, col=fisher.col, pch=19, cex=0.5, main="Bclust", xlab="Easting (m)", ylab="Northing (m)") legend("topleft", fill=attr(fisher.col, "palette"), legend=names(attr(fisher.col, "table")), bty="n") ################################################### ### chunk number 6: DBHClasses ################################################### fixed <- classIntervals(DBH, n=4, style="fixed", fixedBreaks=c(0,12.7,30.48,60,max(DBH)+1)) fixed.col <- findColours(fixed, col.pal) plot(coords, col=fixed.col, pch=19, cex=0.5, main="Forestry tree size classes", xlab="Easting (m)", ylab="Northing (m)") legend("topleft", fill=attr(fixed.col, "palette"), legend=c("sapling","poletimber","sawtimber","large sawtimber"), bty="n") ################################################### ### chunk number 7: MBASurf ################################################### x.res <- 100; y.res <- 100 surf <- mba.surf(cbind(coords, DBH), no.X=x.res, no.Y=y.res, h=5, m=2, extend=FALSE)$xyz.est image.plot(surf, xaxs = "r", yaxs = "r", xlab="Easting (m)", ylab="Northing (m)", col=col.br(25)) ################################################### ### chunk number 8: rglMBASurf eval=FALSE ################################################### ## col <- rbind(0, cbind(matrix(drape.color(surf[[3]], col=col.br(25)), x.res-1, y.res-1), 0)) ## ## surface3d(surf[[1]], surf[[2]], surf[[3]], col=col) ## axes3d(); title3d(main="DBH", xlab="Easting (m)", ylab="Northing (m)", zlab="DBH (cm)") ################################################### ### chunk number 9: perspMBASurf ################################################### drape.plot(surf[[1]], surf[[2]], surf[[3]], col=col.br(150), theta=225, phi=50, border=FALSE, add.legend=FALSE, xlab="Easting (m)", ylab="Northing (m)", zlab="DBH (cm)") image.plot(zlim=range(surf[[3]], na.rm=TRUE), legend.only=TRUE, horizontal=FALSE) ################################################### ### chunk number 10: geoRVariog ################################################### max.dist <- 0.25*max(iDist(coords)) bins <- 50 vario.DBH <- variog(coords=coords, data=DBH, uvec=(seq(0, max.dist, length=bins))) fit.DBH <-variofit(vario.DBH, ini.cov.pars=c(600, 200/-log(0.05)), ##sigma^2 and 1/phi cov.model="exponential", minimisation.function="nls", weights="equal") lm.DBH <- lm(DBH~Species, data=WEF.dat) summary(lm.DBH) DBH.resid <- resid(lm.DBH) vario.DBH.resid <- variog(coords=coords, data=DBH.resid, uvec=(seq(0, max.dist, length=bins))) fit.DBH.resid <-variofit(vario.DBH.resid, ini.cov.pars=c(300, 200/-log(0.05)), cov.model="exponential", minimisation.function="nls", weights="equal") par(mfrow=c(1,2)) plot(vario.DBH, ylim=c(200,1200), main="DBH") lines(fit.DBH) abline(h=fit.DBH$nugget, col="blue")##nugget abline(h=fit.DBH$cov.pars[1]+fit.DBH$nugget, col="green")##sill abline(v=-log(0.05)*fit.DBH$cov.pars[2], col="red3")##effective range plot(vario.DBH.resid, ylim=c(200,500), main="DBH residuals") lines(fit.DBH.resid) abline(h=fit.DBH.resid$nugget, col="blue") abline(h=fit.DBH.resid$cov.pars[1]+fit.DBH.resid$nugget, col="green") abline(v=-log(0.05)*fit.DBH.resid$cov.pars[2], col="red3") ################################################### ### chunk number 11: anisoVario ################################################### ELEV_m <- WEF.dat$ELEV_m sp.dat <- as.data.frame(cbind(DBH, HT, coords, ELEV_m)) coordinates(sp.dat) <- c("East_m", "North_m") vario.DBH <- variogram(DBH~ELEV_m, data=sp.dat, cutoff=max.dist, width=5, alpha=(0:3)*45) fit.DBH <- fit.variogram(vario.DBH, vgm(1000,"Exp",200/-log(0.05),600))##sill, range, nugget print(plot(vario.DBH, fit.DBH)) ################################################### ### chunk number 12: coVario ################################################### g <- gstat(NULL, "DBH", DBH~ELEV_m, data=sp.dat) g <- gstat(g, "HT", HT~ELEV_m, data=sp.dat) g vm <- variogram(g) vm.fit <- fit.lmc(vm, g, vgm(500, "Exp", 100/-log(0.05), 300)) vm.fit print(plot(vm, vm.fit)) cov(cbind(DBH, HT))