Nothing
inst/doc/Vquantile.R
In ModelMap: Modeling and Map Production using Random Forest and Related Stochastic Models
### R code from vignette source 'Vquantile.Rnw'
###################################################
### code chunk number 1: ExSetup set options
###################################################
options(prompt = "R> ")
options(width = 75)
options(continue=" ")
pdf("Vplots.pdf")
###################################################
### code chunk number 2: ExSetup set width
###################################################
options(width=60)
###################################################
### code chunk number 3: ExSetup load package
###################################################
library("ModelMap")
###################################################
### code chunk number 4: ExElevReadGDAL
###################################################
library(raster)
elevfn <- paste(getwd(),"/VModelMapData_dem_ELEVM_250.img",sep="")
mapgrid <- raster(elevfn)
###################################################
### code chunk number 5: ExElev
###################################################
opar <- par(mar=c(4,4,3,6),xpd=NA,mgp=c(3, 2, .3))
col.ramp<-terrain.colors(101)
zlim <- c(1500,maxValue(mapgrid))
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"m",sep="")
legend.label<-paste((7:3)*500,"m")
legend.colors<-col.ramp[c(100,75,50,25,1)]
image( mapgrid,
col = col.ramp,
xlab="", ylab="",
zlim=zlim,
asp=1, bty="n", main="")
legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("Elevation of Study Region",side=3,line=1,cex=1.5)
par(opar)
###################################################
### code chunk number 6: ExElevFig
###################################################
opar <- par(mar=c(4,4,3,6),xpd=NA,mgp=c(3, 2, .3))
col.ramp<-terrain.colors(101)
zlim <- c(1500,maxValue(mapgrid))
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"m",sep="")
legend.label<-paste((7:3)*500,"m")
legend.colors<-col.ramp[c(100,75,50,25,1)]
image( mapgrid,
col = col.ramp,
xlab="", ylab="",
zlim=zlim,
asp=1, bty="n", main="")
legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("Elevation of Study Region",side=3,line=1,cex=1.5)
par(opar)
###################################################
### code chunk number 7: ExSetup Define training and test files
###################################################
qdatafn <- "VModelMapData.csv"
qdata.trainfn <- "VModelMapData_TRAIN.csv"
qdata.testfn <- "VModelMapData_TEST.csv"
###################################################
### code chunk number 8: ExSetup define folder
###################################################
folder <- getwd()
###################################################
### code chunk number 9: ExSetup split training and test
###################################################
get.test( proportion.test=0.2,
qdatafn=qdatafn,
seed=42,
folder=folder,
qdata.trainfn=qdata.trainfn,
qdata.testfn=qdata.testfn)
###################################################
### code chunk number 10: ExSetup Define predictors
###################################################
predList <- c( "ELEV250",
"NLCD01_250",
"EVI2005097",
"NDV2005097",
"NIR2005097",
"RED2005097")
predFactor <- c("NLCD01_250")
###################################################
### code chunk number 11: Ex1 Define Identifier
###################################################
unique.rowname <- "ID"
###################################################
### code chunk number 12: ExSetup update raster LUT
###################################################
rastLUTfn <- "VModelMapData_LUT.csv"
rastLUTfn <- read.table( rastLUTfn,
header=FALSE,
sep=",",
stringsAsFactors=FALSE)
rastLUTfn[,1] <- paste(folder,rastLUTfn[,1],sep="/")
###################################################
### code chunk number 13: Ex1 Define model filenames
###################################################
MODELfn.pinyon <- "VQuantile_QRF_Pinyon"
MODELfn.sage <- "VQuantile_QRF_Sage"
###################################################
### code chunk number 14: Ex1 Define response
###################################################
response.name.pinyon <- "PINYON"
response.name.sage <- "SAGE"
response.type <- "continuous"
###################################################
### code chunk number 15: Ex1 Create Model
###################################################
QRF.pinyon <- model.build( model.type="QRF",
qdata.trainfn=qdata.trainfn,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.pinyon,
predList=predList,
predFactor=predFactor,
response.name=response.name.pinyon,
response.type=response.type,
#importance currently unavailable for QRF
#importance=TRUE,
#quantiles=c(0.1,0.2,0.5,0.8,0.9)
)
QRF.sage <- model.build( model.type="QRF",
qdata.trainfn=qdata.trainfn,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.sage,
predList=predList,
predFactor=predFactor,
response.name=response.name.sage,
response.type=response.type,
#importance currently unavailable for QRF
#importance=TRUE,
#quantiles=c(0.1,0.2,0.5,0.8,0.9)
)
###################################################
### code chunk number 16: Ex1 Model Diagnostics
###################################################
QRF.pinyon.pred <- model.diagnostics( model.obj=QRF.pinyon,
qdata.testfn=qdata.testfn,
folder=folder,
MODELfn=MODELfn.pinyon,
unique.rowname=unique.rowname,
quantiles=c(0.1,0.5,0.9),
# Model Validation Arguments
prediction.type="TEST",
device.type=c("pdf","png"),
cex=1.2)
QRF.sage.pred <- model.diagnostics( model.obj=QRF.sage,
qdata.testfn=qdata.testfn,
folder=folder,
MODELfn=MODELfn.sage,
unique.rowname=unique.rowname,
quantiles=c(0.1,0.5,0.9),
# Model Validation Arguments
prediction.type="TEST",
device.type=c("pdf","png"),
cex=1.2)
###################################################
### code chunk number 17: Ex1InteractionPlots
###################################################
pred.means <- c( ELEV250 = 2300,
NLCD01_250 = 42,
EVI2005097 = 1800,
NDV2005097 = 3100,
NIR2005097 = 2000,
RED2005097 = 1000)
opar <- par(mfrow=c(3,1),mar=c(2,3,0,2),oma=c(0,0,3,0))
for(quantile in c(0.1,0.5,0.9)){
model.interaction.plot( QRF.pinyon$QRF,
x="ELEV250",
y="EVI2005097",
main="",
plot.type="persp",
device.type=c("none"),
MODELfn=MODELfn.pinyon,
folder=paste(folder,"/interaction",sep=""),
quantile=quantile,
pred.means=pred.means,
zlim=c(0,60))
mtext(paste( quantile*100,"% Quantile",sep=""),side=2,line=1,font=2,cex=1)
}
mtext("Pinyon Cover by Quantile",side=3,line=1,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 18: Ex1InteractionPlots
###################################################
pred.means <- c( ELEV250 = 2300,
NLCD01_250 = 42,
EVI2005097 = 1800,
NDV2005097 = 3100,
NIR2005097 = 2000,
RED2005097 = 1000)
opar <- par(mfrow=c(3,1),mar=c(2,3,0,2),oma=c(0,0,3,0))
for(quantile in c(0.1,0.5,0.9)){
model.interaction.plot( QRF.pinyon$QRF,
x="ELEV250",
y="EVI2005097",
main="",
plot.type="persp",
device.type=c("none"),
MODELfn=MODELfn.pinyon,
folder=paste(folder,"/interaction",sep=""),
quantile=quantile,
pred.means=pred.means,
zlim=c(0,60))
mtext(paste( quantile*100,"% Quantile",sep=""),side=2,line=1,font=2,cex=1)
}
mtext("Pinyon Cover by Quantile",side=3,line=1,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 19: Ex1 update raster LUT
###################################################
rastLUTfn <- "VModelMapData_LUT.csv"
rastLUTfn <- read.table( rastLUTfn,
header=FALSE,
sep=",",
stringsAsFactors=FALSE)
rastLUTfn[,1] <- paste(folder,rastLUTfn[,1],sep="/")
###################################################
### code chunk number 20: Ex1 Produce Maps
###################################################
model.mapmake( model.obj=QRF.pinyon,
folder=folder,
MODELfn=MODELfn.pinyon,
rastLUTfn=rastLUTfn,
na.action="na.omit",
# Mapping arguments
map.sd=TRUE,
quantiles=c(0.1,0.5,0.9))
model.mapmake( model.obj=QRF.sage,
folder=folder,
MODELfn=MODELfn.sage,
rastLUTfn=rastLUTfn,
na.action="na.omit",
# Mapping arguments
map.sd=TRUE,
quantiles=c(0.1,0.5,0.9))
###################################################
### code chunk number 21: Ex1 define color sequence
###################################################
l <- seq(100,0,length.out=101)
c <- seq(0,100,length.out=101)
col.ramp <- hcl(h = 120, c = c, l = l)
###################################################
### code chunk number 22: Ex1RFMap
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map_RF.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map_RF.img",sep=""))
zlim <- c(0,60)
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"%",sep="")
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("RF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 23: Ex1RFMapFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map_RF.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map_RF.img",sep=""))
zlim <- c(0,60)
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"%",sep="")
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("RF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 24: Ex1QRFMap
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- brick(paste(MODELfn.pinyon,"_map_QRF.img",sep=""))
mapgrid.sage <- brick(paste(MODELfn.sage,"_map_QRF.img",sep=""))
image( mapgrid.pinyon[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Median (50% quantile) Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 25: Ex1QRFMapFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- brick(paste(MODELfn.pinyon,"_map_QRF.img",sep=""))
mapgrid.sage <- brick(paste(MODELfn.sage,"_map_QRF.img",sep=""))
image( mapgrid.pinyon[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Median (50% quantile) Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 26: Ex1QRFMapLow
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Lower bound (10% quantile)",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 27: Ex1QRFMapLowFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Lower bound (10% quantile)",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 28: Ex1QRFMapUp
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF model - Upper - 90% quantile",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 29: Ex1QRFMapUpFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF model - Upper - 90% quantile",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 30: Ex1 Define model filenames
###################################################
MODELfn.pinyon <- "VQuantile_CF_Pinyon"
MODELfn.sage <- "VQuantile_CF_Sage"
###################################################
### code chunk number 31: Ex1 Define response
###################################################
response.name.pinyon <- "PINYON"
response.name.sage <- "SAGE"
response.type <- "continuous"
###################################################
### code chunk number 32: Ex1 Create Model
###################################################
qdata<-read.csv(qdata.trainfn)
IS.NUM<-sapply(qdata,is.numeric)
qdata[,IS.NUM]<-sapply(qdata[,IS.NUM],as.numeric)
CF.pinyon <- model.build( model.type="CF",
qdata.trainfn=qdata,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.pinyon,
predList=predList,
predFactor=predFactor,
response.name=response.name.pinyon,
response.type=response.type)
CF.sage <- model.build( model.type="CF",
qdata.trainfn=qdata,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.sage,
predList=predList,
predFactor=predFactor,
response.name=response.name.sage,
response.type=response.type)
###################################################
### code chunk number 33: Ex2CompImpCF_RF
###################################################
opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))
Imp.pinyon<-model.importance.plot( model.obj.1=CF.pinyon,
model.obj.2=QRF.pinyon$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Pinyon Percent Cover",
device.type="none",
cex=0.9)
Imp.sage<-model.importance.plot( model.obj.1=CF.sage,
model.obj.2=QRF.sage$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Sage Percent Cover",
device.type="none",
cex=0.9)
mtext("CF versus RF Variable Importance",side=3,line=0,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 34: Ex2CompImpFigCF_RF
###################################################
opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))
Imp.pinyon<-model.importance.plot( model.obj.1=CF.pinyon,
model.obj.2=QRF.pinyon$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Pinyon Percent Cover",
device.type="none",
cex=0.9)
Imp.sage<-model.importance.plot( model.obj.1=CF.sage,
model.obj.2=QRF.sage$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Sage Percent Cover",
device.type="none",
cex=0.9)
mtext("CF versus RF Variable Importance",side=3,line=0,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 35: Ex1 Produce Maps
###################################################
model.mapmake( model.obj=CF.pinyon,
folder=folder,
MODELfn=MODELfn.pinyon,
rastLUTfn=rastLUTfn,
na.action="na.omit"
)
model.mapmake( model.obj=CF.sage,
folder=folder,
MODELfn=MODELfn.sage,
rastLUTfn=rastLUTfn,
na.action="na.omit"
)
###################################################
### code chunk number 36: Ex2CFMap
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map.img",sep=""))
zlim <- c(0,60)
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("CF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 37: Ex2CFMapFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map.img",sep=""))
zlim <- c(0,60)
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("CF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 38: Remove Rplots pdf
###################################################
dev.off()
file.remove("Vplots.pdf")
Try the ModelMap package in your browser
Any scripts or data that you put into this service are public.
ModelMap documentation built on April 25, 2023, 1:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### R code from vignette source 'Vquantile.Rnw'
###################################################
### code chunk number 1: ExSetup set options
###################################################
options(prompt = "R> ")
options(width = 75)
options(continue=" ")
pdf("Vplots.pdf")
###################################################
### code chunk number 2: ExSetup set width
###################################################
options(width=60)
###################################################
### code chunk number 3: ExSetup load package
###################################################
library("ModelMap")
###################################################
### code chunk number 4: ExElevReadGDAL
###################################################
library(raster)
elevfn <- paste(getwd(),"/VModelMapData_dem_ELEVM_250.img",sep="")
mapgrid <- raster(elevfn)
###################################################
### code chunk number 5: ExElev
###################################################
opar <- par(mar=c(4,4,3,6),xpd=NA,mgp=c(3, 2, .3))
col.ramp<-terrain.colors(101)
zlim <- c(1500,maxValue(mapgrid))
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"m",sep="")
legend.label<-paste((7:3)*500,"m")
legend.colors<-col.ramp[c(100,75,50,25,1)]
image( mapgrid,
col = col.ramp,
xlab="", ylab="",
zlim=zlim,
asp=1, bty="n", main="")
legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("Elevation of Study Region",side=3,line=1,cex=1.5)
par(opar)
###################################################
### code chunk number 6: ExElevFig
###################################################
opar <- par(mar=c(4,4,3,6),xpd=NA,mgp=c(3, 2, .3))
col.ramp<-terrain.colors(101)
zlim <- c(1500,maxValue(mapgrid))
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"m",sep="")
legend.label<-paste((7:3)*500,"m")
legend.colors<-col.ramp[c(100,75,50,25,1)]
image( mapgrid,
col = col.ramp,
xlab="", ylab="",
zlim=zlim,
asp=1, bty="n", main="")
legend( x=xmax(mapgrid),y=ymax(mapgrid),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("Elevation of Study Region",side=3,line=1,cex=1.5)
par(opar)
###################################################
### code chunk number 7: ExSetup Define training and test files
###################################################
qdatafn <- "VModelMapData.csv"
qdata.trainfn <- "VModelMapData_TRAIN.csv"
qdata.testfn <- "VModelMapData_TEST.csv"
###################################################
### code chunk number 8: ExSetup define folder
###################################################
folder <- getwd()
###################################################
### code chunk number 9: ExSetup split training and test
###################################################
get.test( proportion.test=0.2,
qdatafn=qdatafn,
seed=42,
folder=folder,
qdata.trainfn=qdata.trainfn,
qdata.testfn=qdata.testfn)
###################################################
### code chunk number 10: ExSetup Define predictors
###################################################
predList <- c( "ELEV250",
"NLCD01_250",
"EVI2005097",
"NDV2005097",
"NIR2005097",
"RED2005097")
predFactor <- c("NLCD01_250")
###################################################
### code chunk number 11: Ex1 Define Identifier
###################################################
unique.rowname <- "ID"
###################################################
### code chunk number 12: ExSetup update raster LUT
###################################################
rastLUTfn <- "VModelMapData_LUT.csv"
rastLUTfn <- read.table( rastLUTfn,
header=FALSE,
sep=",",
stringsAsFactors=FALSE)
rastLUTfn[,1] <- paste(folder,rastLUTfn[,1],sep="/")
###################################################
### code chunk number 13: Ex1 Define model filenames
###################################################
MODELfn.pinyon <- "VQuantile_QRF_Pinyon"
MODELfn.sage <- "VQuantile_QRF_Sage"
###################################################
### code chunk number 14: Ex1 Define response
###################################################
response.name.pinyon <- "PINYON"
response.name.sage <- "SAGE"
response.type <- "continuous"
###################################################
### code chunk number 15: Ex1 Create Model
###################################################
QRF.pinyon <- model.build( model.type="QRF",
qdata.trainfn=qdata.trainfn,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.pinyon,
predList=predList,
predFactor=predFactor,
response.name=response.name.pinyon,
response.type=response.type,
#importance currently unavailable for QRF
#importance=TRUE,
#quantiles=c(0.1,0.2,0.5,0.8,0.9)
)
QRF.sage <- model.build( model.type="QRF",
qdata.trainfn=qdata.trainfn,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.sage,
predList=predList,
predFactor=predFactor,
response.name=response.name.sage,
response.type=response.type,
#importance currently unavailable for QRF
#importance=TRUE,
#quantiles=c(0.1,0.2,0.5,0.8,0.9)
)
###################################################
### code chunk number 16: Ex1 Model Diagnostics
###################################################
QRF.pinyon.pred <- model.diagnostics( model.obj=QRF.pinyon,
qdata.testfn=qdata.testfn,
folder=folder,
MODELfn=MODELfn.pinyon,
unique.rowname=unique.rowname,
quantiles=c(0.1,0.5,0.9),
# Model Validation Arguments
prediction.type="TEST",
device.type=c("pdf","png"),
cex=1.2)
QRF.sage.pred <- model.diagnostics( model.obj=QRF.sage,
qdata.testfn=qdata.testfn,
folder=folder,
MODELfn=MODELfn.sage,
unique.rowname=unique.rowname,
quantiles=c(0.1,0.5,0.9),
# Model Validation Arguments
prediction.type="TEST",
device.type=c("pdf","png"),
cex=1.2)
###################################################
### code chunk number 17: Ex1InteractionPlots
###################################################
pred.means <- c( ELEV250 = 2300,
NLCD01_250 = 42,
EVI2005097 = 1800,
NDV2005097 = 3100,
NIR2005097 = 2000,
RED2005097 = 1000)
opar <- par(mfrow=c(3,1),mar=c(2,3,0,2),oma=c(0,0,3,0))
for(quantile in c(0.1,0.5,0.9)){
model.interaction.plot( QRF.pinyon$QRF,
x="ELEV250",
y="EVI2005097",
main="",
plot.type="persp",
device.type=c("none"),
MODELfn=MODELfn.pinyon,
folder=paste(folder,"/interaction",sep=""),
quantile=quantile,
pred.means=pred.means,
zlim=c(0,60))
mtext(paste( quantile*100,"% Quantile",sep=""),side=2,line=1,font=2,cex=1)
}
mtext("Pinyon Cover by Quantile",side=3,line=1,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 18: Ex1InteractionPlots
###################################################
pred.means <- c( ELEV250 = 2300,
NLCD01_250 = 42,
EVI2005097 = 1800,
NDV2005097 = 3100,
NIR2005097 = 2000,
RED2005097 = 1000)
opar <- par(mfrow=c(3,1),mar=c(2,3,0,2),oma=c(0,0,3,0))
for(quantile in c(0.1,0.5,0.9)){
model.interaction.plot( QRF.pinyon$QRF,
x="ELEV250",
y="EVI2005097",
main="",
plot.type="persp",
device.type=c("none"),
MODELfn=MODELfn.pinyon,
folder=paste(folder,"/interaction",sep=""),
quantile=quantile,
pred.means=pred.means,
zlim=c(0,60))
mtext(paste( quantile*100,"% Quantile",sep=""),side=2,line=1,font=2,cex=1)
}
mtext("Pinyon Cover by Quantile",side=3,line=1,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 19: Ex1 update raster LUT
###################################################
rastLUTfn <- "VModelMapData_LUT.csv"
rastLUTfn <- read.table( rastLUTfn,
header=FALSE,
sep=",",
stringsAsFactors=FALSE)
rastLUTfn[,1] <- paste(folder,rastLUTfn[,1],sep="/")
###################################################
### code chunk number 20: Ex1 Produce Maps
###################################################
model.mapmake( model.obj=QRF.pinyon,
folder=folder,
MODELfn=MODELfn.pinyon,
rastLUTfn=rastLUTfn,
na.action="na.omit",
# Mapping arguments
map.sd=TRUE,
quantiles=c(0.1,0.5,0.9))
model.mapmake( model.obj=QRF.sage,
folder=folder,
MODELfn=MODELfn.sage,
rastLUTfn=rastLUTfn,
na.action="na.omit",
# Mapping arguments
map.sd=TRUE,
quantiles=c(0.1,0.5,0.9))
###################################################
### code chunk number 21: Ex1 define color sequence
###################################################
l <- seq(100,0,length.out=101)
c <- seq(0,100,length.out=101)
col.ramp <- hcl(h = 120, c = c, l = l)
###################################################
### code chunk number 22: Ex1RFMap
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map_RF.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map_RF.img",sep=""))
zlim <- c(0,60)
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"%",sep="")
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("RF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 23: Ex1RFMapFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map_RF.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map_RF.img",sep=""))
zlim <- c(0,60)
legend.label<-rev(pretty(zlim,n=5))
legend.colors<-col.ramp[trunc((legend.label/max(legend.label))*100)+1]
legend.label<-paste(legend.label,"%",sep="")
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("RF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 24: Ex1QRFMap
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- brick(paste(MODELfn.pinyon,"_map_QRF.img",sep=""))
mapgrid.sage <- brick(paste(MODELfn.sage,"_map_QRF.img",sep=""))
image( mapgrid.pinyon[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Median (50% quantile) Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 25: Ex1QRFMapFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- brick(paste(MODELfn.pinyon,"_map_QRF.img",sep=""))
mapgrid.sage <- brick(paste(MODELfn.sage,"_map_QRF.img",sep=""))
image( mapgrid.pinyon[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[2]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Median (50% quantile) Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 26: Ex1QRFMapLow
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Lower bound (10% quantile)",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 27: Ex1QRFMapLowFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[1]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF - Lower bound (10% quantile)",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 28: Ex1QRFMapUp
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF model - Upper - 90% quantile",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 29: Ex1QRFMapUpFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
image( mapgrid.pinyon[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage[[3]],
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("QRF model - Upper - 90% quantile",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 30: Ex1 Define model filenames
###################################################
MODELfn.pinyon <- "VQuantile_CF_Pinyon"
MODELfn.sage <- "VQuantile_CF_Sage"
###################################################
### code chunk number 31: Ex1 Define response
###################################################
response.name.pinyon <- "PINYON"
response.name.sage <- "SAGE"
response.type <- "continuous"
###################################################
### code chunk number 32: Ex1 Create Model
###################################################
qdata<-read.csv(qdata.trainfn)
IS.NUM<-sapply(qdata,is.numeric)
qdata[,IS.NUM]<-sapply(qdata[,IS.NUM],as.numeric)
CF.pinyon <- model.build( model.type="CF",
qdata.trainfn=qdata,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.pinyon,
predList=predList,
predFactor=predFactor,
response.name=response.name.pinyon,
response.type=response.type)
CF.sage <- model.build( model.type="CF",
qdata.trainfn=qdata,
folder=folder,
unique.rowname=unique.rowname,
MODELfn=MODELfn.sage,
predList=predList,
predFactor=predFactor,
response.name=response.name.sage,
response.type=response.type)
###################################################
### code chunk number 33: Ex2CompImpCF_RF
###################################################
opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))
Imp.pinyon<-model.importance.plot( model.obj.1=CF.pinyon,
model.obj.2=QRF.pinyon$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Pinyon Percent Cover",
device.type="none",
cex=0.9)
Imp.sage<-model.importance.plot( model.obj.1=CF.sage,
model.obj.2=QRF.sage$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Sage Percent Cover",
device.type="none",
cex=0.9)
mtext("CF versus RF Variable Importance",side=3,line=0,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 34: Ex2CompImpFigCF_RF
###################################################
opar <- par(mfrow=c(2,1),mar=c(3,3,3,3),oma=c(0,0,3,0))
Imp.pinyon<-model.importance.plot( model.obj.1=CF.pinyon,
model.obj.2=QRF.pinyon$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Pinyon Percent Cover",
device.type="none",
cex=0.9)
Imp.sage<-model.importance.plot( model.obj.1=CF.sage,
model.obj.2=QRF.sage$RF,
model.name.1="unconditional (CF)",
model.name.2="unconditional (RF)",
imp.type.1=1,
imp.type.2=2,
cf.conditional.1=FALSE,
sort.by="predList",
predList=predList,
scale.by="sum",
main="Sage Percent Cover",
device.type="none",
cex=0.9)
mtext("CF versus RF Variable Importance",side=3,line=0,cex=1.8,outer=TRUE)
par(opar)
###################################################
### code chunk number 35: Ex1 Produce Maps
###################################################
model.mapmake( model.obj=CF.pinyon,
folder=folder,
MODELfn=MODELfn.pinyon,
rastLUTfn=rastLUTfn,
na.action="na.omit"
)
model.mapmake( model.obj=CF.sage,
folder=folder,
MODELfn=MODELfn.sage,
rastLUTfn=rastLUTfn,
na.action="na.omit"
)
###################################################
### code chunk number 36: Ex2CFMap
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map.img",sep=""))
zlim <- c(0,60)
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("CF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 37: Ex2CFMapFig
###################################################
opar <- par(mfrow=c(1,2),mar=c(3,3,2,1),oma=c(0,0,3,4),xpd=NA)
mapgrid.pinyon <- raster(paste(MODELfn.pinyon,"_map.img",sep=""))
mapgrid.sage <- raster(paste(MODELfn.sage,"_map.img",sep=""))
zlim <- c(0,60)
image( mapgrid.pinyon,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.pinyon,side=3,line=1,cex=1.2)
image( mapgrid.sage,
col=col.ramp,
xlab="",ylab="",xaxt="n",yaxt="n",
zlim=zlim,
asp=1,bty="n",main="")
mtext(response.name.sage,side=3,line=1,cex=1.2)
legend( x=xmax(mapgrid.sage),y=ymax(mapgrid.sage),
legend=legend.label,
fill=legend.colors,
bty="n",
cex=1.2)
mtext("CF - Mean Percent Cover",side=3,line=1,cex=1.5,outer=T)
par(opar)
###################################################
### code chunk number 38: Remove Rplots pdf
###################################################
dev.off()
file.remove("Vplots.pdf")
Try the ModelMap package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.