library(spBayes)
library(classInt)
library(RColorBrewer)
library(MBA)
library(fields)
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.6-0 (2020-12-14) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: viridis
## Loading required package: viridisLite
## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code
library(rgl)
data("WEF.dat")
head(WEF.dat)
## Tree_ID Species DBH_cm Live_dead Stump North_m East_m ELEV_m Tree_height_m
## 1 5 SF 63.3 D 79.03 221.85 1112.15 NA
## 2 18 SF 38.0 D 100.70 225.05 1113.03 NA
## 3 20 SF 16.4 D 95.41 223.93 1112.76 NA
## 4 26 WH 12.6 D 108.01 231.90 1113.46 NA
## 5 30 SF 21.5 D 98.43 233.42 1112.80 NA
## 6 39 SF 42.4 D Stump 88.15 239.38 1113.20 NA
## Crown_depth_m Mean_crown_radius_m_from_center
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
## Crown_radius_north_m_from_center Crown_radius_east_m_from_center
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
## Crown_radius_west_m_from_center Crown_radius_south_m_from_center
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
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(m)")
col.br <- colorRampPalette(c("blue","cyan","yellow","red"))
col.pal <- col.br(5)
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=1,cex=0.5,
main="Forest tree size classes",
xlab="Easting(m)",ylab="Northing(m)")
legend("topleft",fill = attr(fixed.col,"palette"),
legend = c("sapling","poletimeber","sawtimber","large sawtimber"),bty = "n")
y.res<- 100
x.res<- 100
surf <- mba.surf(cbind(coords,DBH),no.X = x.res,no.Y = y.res,
h=5, m = 2, extend = FALSE)
surf <- surf$xyz.est
image.plot(surf,xaxs = "r",yaxs = "r",xlab="Easting(m)",ylab="Northing(m)",col = col.br(25))
contour(surf,add=T)
#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)")
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 )
