summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--figures.R388
1 files changed, 388 insertions, 0 deletions
diff --git a/figures.R b/figures.R
new file mode 100644
index 0000000..1b9c0c5
--- /dev/null
+++ b/figures.R
@@ -0,0 +1,388 @@
+
+#Data sources
+source('data-sources.R')
+
+#Load required packages
+if(!'png'%in%installed.packages()[,'Package']){
+ install.packages('png',repos="https://cran.mtu.edu")}
+
+#Color palette
+cols.mu = c('black','firebrick1','dodgerblue')
+cols.hl <- c('black','firebrick1','orange','dodgerblue','forestgreen')
+
+# BACI figure
+
+#Simulate data
+set.seed(2016)
+control <- rnorm(10,0.3,0.05)
+treatment <- c(rnorm(3,0.3,0.05),rnorm(7,1,0.1))
+
+pdf("figures/baci.pdf",width=5.7,height=4, pointsize=10)
+
+par(mgp=c(3,1,0),mar = c(4.5,4.5,1,1) + 0.1,oma=c(0,0,0,0))
+structure <- c(1:10)
+plot(structure,treatment,type='l',col=cols.mu[2],lwd=2,ylim=c(0,1.5),
+ xlab='Year',ylab='Density', cex.lab=1.7)
+
+lines(structure,control,col=cols.mu[1],lwd=2)
+abline(v=3.5,lty=2)
+
+text(3,1,'Treatment Applied',srt=90,cex=1.7)
+text(mean(c(0.5,3.5)),0.1,'Before',font=2,cex=1.7)
+text(mean(c(3.5,10.5)),0.1,'After',font=2,cex=1.7)
+
+arrows(7,0.7,7,control[7]+0.05,length=0.1,lwd=2)
+arrows(7,0.7,7,treatment[7]-0.05,length=0.1,lwd=2)
+
+arrows(2,0.33,2,control[2]+0.01,length=0.05,lwd=2)
+arrows(2,0.33,2,treatment[2]-0.01,length=0.05,lwd=2)
+
+legend('topright',col=cols.mu[2:1],lwd=2,legend=c('Treatment','Control'),cex=1.5)
+
+dev.off()
+
+#######################################
+
+#Species response, harvest scale
+source("https://raw.githubusercontent.com/kenkellner/hee-bird-abundance/master/figures/figure_resp_summary.R")
+
+pdf("figures/hl_resp.pdf",width=6.5,height=4,pointsize=10)
+
+resp.summary.figure('hl',c(-8.5,20.5),cols.hl[c(2,2,3,3)],text.cex=0.65,
+ highlight.shrubland=TRUE,lab.cex=1.3)
+
+dev.off()
+
+##################################################
+#Individual Species, Harvest Scale
+source('https://raw.githubusercontent.com/kenkellner/hee-bird-abundance/master/figures/figure_birddens.R')
+
+pdf("figures/focal_hl.pdf",width=6.5,height=4,pointsize=9)
+
+par(fig=c(0,0.54,0.43,1),new=FALSE)
+bird.figure('BAWW','hl',color=cols.hl,title=
+ 'Black-and-white Warbler',legendpos='none',legcex=1,
+ extralabels=F,xlab=NULL,se=F,ylab="")
+par(fig=c(0.46,1,0.43,1),new=TRUE)
+bird.figure('HOWA','hl',color=cols.hl,title=
+ 'Hooded Warbler',legendpos='none',legcex=1,
+ extralabels=F,xlab=NULL,ylab="",se=F)
+par(fig=c(0,0.54,0,0.57),new=TRUE)
+bird.figure('CERW','hl',color=cols.hl,title=
+ 'Cerulean Warbler',legendpos='topright',legcex=0.8,
+ legend.lwd=1, extralabels=F,xlab="",se=F,ylab="")
+par(fig=c(0.46,1,0,0.57),new=TRUE)
+bird.figure('WEWA','hl',color=cols.hl,title=
+ 'Worm-eating Warbler',legendpos='none',legcex=1,
+ extralabels=F,xlab="",ylab="",se=F)
+mtext("Density (males/ha)",2,outer=T,line=-2,cex=2)
+dev.off()
+
+#################################
+
+# Harvest-level Bird Richness figure
+source("https://raw.githubusercontent.com/kenkellner/hee-bird-abundance/master/figures/figure_richness.R")
+
+pdf("figures/hl_rich.pdf",width=6,height=4,pointsize=10)
+richness.figure(scale='hl',figmin=16.5,figmax=32,cols=cols.hl,
+ pchs=c(17,15,16,18,21),ltys=rep(1,5),xlab="",lab.cex=1.7,leg.cex=1.3)
+dev.off()
+
+##
+#####################################################
+#Rubus/SCTA Figure
+
+#Read in SCTA data
+birds$Date <- as.Date(as.character(birds$Date), "%m/%d/%y")
+birds.scta <- birds[birds$Species=='Scarlet Tanager',]
+
+#Organize by date/session
+scta.sum <- aggregate(birds.scta$Age,
+ list(Date=birds.scta$Session, Location=birds.scta$Location),length,
+ drop=FALSE)
+scta.sum$Date <- aggregate(birds.scta$Date,list(Session=birds.scta$Session,
+ Location=birds.scta$Location),mean,drop=FALSE)$x
+scta.sum <- scta.sum[order(scta.sum$Date),]
+
+#Read in veg data
+veg$Date <- as.Date(as.character(veg$Date), "%m/%d/%Y")
+veg.sort <- veg[order(veg$Date),]
+
+pdf("figures/clearcut_fruit.pdf",width=6.5,height=4,pointsize=10)
+
+#Set up plot
+par(mgp=c(3,1,0),mar = c(4.5,4.5,1,4.5) + 0.1,oma=c(0,0,0,0))
+plot(veg.sort$Date,veg.sort$Rubus,type='l',lwd=3,col="black",
+ xlab='Date',ylab='Rubus Index',cex.lab=1.7)
+par(new=TRUE)
+plot(scta.sum$Date,scta.sum$x,type='l',lwd=3,col="firebrick1",xaxt='n',yaxt='n',
+ xlab='',ylab="")
+axis(4)
+mtext("Scarlet Tanager Captures",side=4,line=3,cex=1.7)
+
+legend('topleft',legend=c('Rubus Index','Scarlet Tanager Captures'),
+ col=cols.mu[1:2],lwd=2)
+
+require(png)
+rub <- readPNG('images/rubus.png')
+tan <- readPNG('images/scarlet_tanager.png')
+rasterImage(rub, scta.sum$Date[1],5,scta.sum$Date[1]+25,10)
+rasterImage(tan,scta.sum$Date[22]-1,10,scta.sum$Date[22]+15,13)
+
+dev.off()
+
+#########################
+
+## Small mammal microhabitat
+
+cols.sm <- c("black","dodgerblue","orange","firebrick1")
+
+clearmean = matrix(NA,nrow=5,ncol=5)
+patchmean = matrix(NA,nrow=5,ncol=5)
+sheltermean = matrix(NA,nrow=5,ncol=5)
+conmean = matrix(NA,nrow=5,ncol=5)
+for (i in 2:5){
+for (j in 1:5){
+conmean[i,j] <- mean(micromn[micromn$year==j&micromn$treat==1,i])
+clearmean[i,j] <- mean(micromn[micromn$year==j&micromn$treat==2,i])
+patchmean[i,j] <- mean(micromn[micromn$year==j&micromn$treat%in%c(3,4),i])
+sheltermean[i,j] <- mean(micromn[micromn$year==j&micromn$treat==5,i])
+}
+}
+conmean[1,] <- c(8.98,16.54,5.75,0.44,19.30)
+sheltermean[1,] <- c(5.46,17.36,3.04,2.06,15.54)
+patchmean[1,] <- c(14.42,20.62,0,0.31,20.78)
+clearmean[1,] <- c(2.68,7.53,4.37,0,24.88)
+
+alldata <- array(NA,dim=c(5,5,4))
+alldata[,,1] <- conmean
+alldata[,,2] <- sheltermean
+alldata[,,3] <- patchmean
+alldata[,,4] <- clearmean
+
+#########################
+
+#Microhabitat subset figure
+pdf("figures/smam_microhab_subset.pdf",width=7,height=4,pointsize=10)
+
+par(mfrow=c(2,3), oma = c(2,2,1,1), mar=c(2,2.5,2.5,1))
+
+ymax <- c(25,25,25,2,4)
+trt <- c('Control','Shelterwood','Patch Cut','Clearcut')
+
+for (i in 1:2){
+
+if(i==1){vars<-c("",'Herb Cover (%)','Woody Cover (%)','CWD (m)',"")
+}else{vars <- rep("",3)}
+
+for (j in 2:4){
+
+barplot(alldata[j,1:5,i],col=cols.sm[i],ylim=c(0,ymax[j]),cex.lab=1.5,main=vars[j],cex.main=1.7)
+abline(v=2.52,lty=2)
+
+if(j==2){mtext(trt[i],side=2,line=2.5,cex=1.7)}
+
+}
+
+}
+
+mtext("Year",side=1,outer=T,line=0.5,cex=1.7)
+#mtext("Harvest Type",side=2,outer=T,line=0.5)
+
+dev.off()
+
+##############################
+
+## SDE Figure
+
+#Read in necessary function
+source("https://raw.githubusercontent.com/kenkellner/hee-acorn-dispersal/master/function_SDE_bootstrap.R")
+
+#Function to make each yearly plot
+gen.sde.fig = function(year,xaxis=F,yaxt=TRUE,legend=F){
+
+ if(year=='y12'){yindex=c(2,11,1,9); yindexw=c(3,11,2,8)}
+ if(year=='y13'){yindex=c(3,12,5,8); yindexw=c(5,10,4,9)}
+ if(year=='y14'){yindex=c(6,10,4,7); yindexw=c(6,12,1,7)}
+
+ #Run bootstrap sims to calculated SDE for black oak w/in year
+ out.bo.con.sq <- sim.prob(1000,tag.data.bo[tag.data.bo$year==year,],
+ out.bo.yr.new.100,'squirrel','control',yindex[1])
+ out.bo.sh.sq <- sim.prob(1000,tag.data.bo[tag.data.bo$year==year,],
+ out.bo.yr.new.100,'squirrel','shelter',yindex[2])
+ p1 <- mean(abs(out.bo.con.sq$sde)<abs(out.bo.sh.sq$sde)) #Check for sig diff
+
+ out.bo.con.de <- sim.prob(1000,tag.data.bo[tag.data.bo$year==year,],
+ out.bo.yr.new.100,'deer','control',yindex[3])
+ out.bo.sh.de <- sim.prob(1000,tag.data.bo[tag.data.bo$year==year,],
+ out.bo.yr.new.100,'deer','shelter',yindex[4])
+ p3 <- mean(abs(out.bo.con.de$sde)<abs(out.bo.sh.de$sde)) #Check for sig diff
+
+ #Also run bootstrap for white oak if year != 2011 (which had no white oak data)
+ out.wo.con.sq <- sim.prob(1000,tag.data.wo,out.wo.yr.new.100,'squirrel','control',yindexw[1])
+ out.wo.sh.sq <- sim.prob(1000,tag.data.wo,out.wo.yr.new.100,'squirrel','shelter',yindexw[2])
+ p2 <- mean(abs(out.wo.con.sq$sde)<abs(out.wo.sh.sq$sde))
+ out.wo.con.de <- sim.prob(1000,tag.data.wo,out.wo.yr.new.100,'deer','control',yindexw[3])
+ out.wo.sh.de <- sim.prob(1000,tag.data.wo,out.wo.yr.new.100,'deer','shelter',yindexw[4])
+ p4 <- mean(abs(out.wo.con.de$sde)<abs(out.wo.sh.de$sde))
+
+ fig.means <- c(mean(out.bo.con.sq$sde),mean(out.bo.sh.sq$sde),
+ mean(out.wo.con.sq$sde),mean(out.wo.sh.sq$sde),
+ mean(out.bo.con.de$sde),mean(out.bo.sh.de$sde),
+ mean(out.wo.con.de$sde),mean(out.wo.sh.de$sde))
+
+ fig.sd <- c(sd(out.bo.con.sq$sde),sd(out.bo.sh.sq$sde),
+ sd(out.wo.con.sq$sde),sd(out.wo.sh.sq$sde),
+ sd(out.bo.con.de$sde),sd(out.bo.sh.de$sde),
+ sd(out.wo.con.de$sde),sd(out.wo.sh.de$sde))
+
+ if(year=='y13'){fig.means[2] <- 0; fig.sd[2] <- 0}
+ if(year=='y14'){fig.means[c(2,3,4,6,7,8)] <- 0; fig.sd[c(2,3,4,6,7,8)] <- 0;}
+
+ #Set up basic plot
+ structure <- c(1,2,3.5,4.5,6.5,7.5,9,10)
+
+ plot(structure,fig.means,ylim=c(0,0.7),xaxt='n',yaxt='n',ylab="",xlab="",
+ pch=rep(c(21,21),4),cex=2,xlim=c(0.5,10.5),bg=rep(c('black','dodgerblue'),4))
+ if(yaxt){
+ axis(2,at=c(0,0.2,0.4,0.6),c('0','0.2','0.4','0.6'),cex.axis=1.7)
+ }
+ abline(v=5.5)
+ if(xaxis){axis(1,at=c(1.5,4,7,9.5),labels=rep(c('BO','WO'),2),tick=F,cex.axis=1.7)}
+
+ #Draw 95% CIs
+ width = 0.1
+ for (i in 1:8){
+ segments(structure[i],max(fig.means[i]-1.96*fig.sd[i],0),
+ structure[i],fig.means[i]+1.96*fig.sd[i])
+ segments(structure[i]-width,fig.means[i]+1.96*fig.sd[i],
+ structure[i]+width,fig.means[i]+1.96*fig.sd[i])
+ segments(structure[i]-width,max(fig.means[i]-1.96*fig.sd[i],0),
+ structure[i]+width,max(fig.means[i]-1.96*fig.sd[i],0))
+ }
+
+ #Draw contrast lines
+ for (i in 1:4){
+ index1 <- i*2 - 1
+ index2 <- i*2
+ segments(structure[index1],fig.means[index1],
+ structure[index2],fig.means[index2],lty=3)
+ }
+
+ points(structure,fig.means,pch=rep(c(21,21),4),cex=3,bg=rep(c('black','dodgerblue'),4))
+
+ #Identify plot sections
+ text(2.75,0.68,'Mouse',cex=2)
+ text(8.25,0.68,'Squirrel',cex=2)
+
+ #Indicate significant differences
+ if(p1<0.05){text(1.5,0.6,'*',cex=2)}
+ if(p2<0.05){text(4,0.6,'*',cex=2)}
+ if(p3<0.05){text(7.5,0.6,'*',cex=2)}
+ if(p4<0.05){text(9.5,0.6,'*',cex=2)}
+
+}
+
+pdf('figures/sde.pdf',width=6.3,height=4,pointsize=10)
+par(mar = c(4.5,2.5,2,0.5) + 0.1,oma=c(1.5,2.5,0,0))
+par(mfrow=c(1,3))
+
+gen.sde.fig('y12',yaxt=TRUE,xaxis=T)
+mtext('2011',3,0,cex=1.5)
+gen.sde.fig('y13',yaxt=FALSE,xaxis=T)
+mtext('2012',3,0,cex=1.5)
+gen.sde.fig('y14',xaxis=T,yaxt=FALSE)
+mtext('2013',3,0,cex=1.5)
+legend(1.7,0.5,legend=c('Control','Shelterwood'),pch=c(21,21),pt.bg=c('black','dodgerblue'),
+ bg='white',cex=1.7,pt.cex=3)
+mtext(expression("Seed Dispersal Effectiveness"),2,0.5,outer=T,cex=1.5)
+mtext(expression("Acorn Species"),1,-0.5,outer=T,cex=1.5)
+
+dev.off()
+
+###############################
+
+## Predation Case Study figure from SOEL
+
+source("https://raw.githubusercontent.com/kenkellner/SOEL/master/utility_functions.R")
+#Generate dataset from SOEL output files
+datalist <- list(avg=weevil.dispersal.average,trt=weevil.dispersal.treateff,yrly=weevil.dispersal.yearlyeff,
+ treat.yrly=weevil.dispersal.treatyearlyeff)
+datalist <- add.newseedlings(datalist,30,38)
+datalist <- add.seedorigin(datalist)
+datalist <- add.acornsum(datalist,30,37)
+datalist <- add.pctgermmean(datalist,30,37)
+
+structure <- c(1,2)
+
+
+cols <- rev(c(rgb(red=244,green=125,blue=66, maxColorValue=255),
+ rgb(red=141,green=213,blue=18, maxColorValue=255)))
+
+pdf("figures/soel_predation.pdf",width=6.7,height=4, pointsize=10)
+
+############################################################
+par(mar = c(2.5,4.5,1,2) + 0.1)
+par(fig=c(0,0.51,0,1),mgp=c(2.5,1,0),new=FALSE)
+h = gen.dataset(datalist,'seedlingsum')
+h$scenario = factor(h$scenario,c('avg','trt','yrly','treat.yrly'))
+
+mns <- aggregate(x=h$seedlingsum/1000,by=list(h$scenario,h$harvest),FUN=mean)[,3][c(1,2,5,6)]
+sds <- aggregate(x=h$seedlingsum/1000,by=list(h$scenario,h$harvest),FUN=sd)[,3][c(1,2,5,6)]
+
+mns <- mns[3:4]
+sds <- sds[3:4]
+
+uplim <- mns+sds
+lowlim <- mns-sds
+
+plot(1,xlim=c(0.5,2.5),ylim=c(.95*min(lowlim),1.1*max(uplim)),xaxt='n',xlab="",
+ ylab=expression("New Seedlings x 1000"),cex.lab=1.5,
+ main='')
+
+wd=0.4
+for(i in 1:2){
+ segments(structure[i],lowlim[i],structure[i],uplim[i])
+ segments(structure[i]+wd/2,lowlim[i],structure[i]-wd/2,lowlim[i])
+ segments(structure[i]+wd/2,uplim[i],structure[i]-wd/2,uplim[i])
+}
+
+points(structure,mns,cex=4,bg=cols,pch=21)
+
+text(structure,(uplim+0.3),c('A','B'))
+
+legend('topleft',pch=21,pt.cex=2,pt.bg=cols,bg='white',title="Simulation Scenario",
+ legend=c('Harvest doesn\'t affect dispersal','Harvest affects dispersal'),cex=1)
+
+#############################################################################################
+
+#Plot number of saplings x scenario
+par(fig=c(0.49,1,0,1),mgp=c(2.5,1,0),new=TRUE)
+h = gen.dataset(datalist,'seedorigin',37)
+h$scenario = factor(h$scenario,c('avg','trt','yrly','treat.yrly'))
+
+mns <- aggregate(x=h$seedorigin,by=list(h$scenario,h$harvest),FUN=mean)[,3][c(1,2,5,6)]
+sds <- aggregate(x=h$seedorigin,by=list(h$scenario,h$harvest),FUN=sd)[,3][c(1,2,5,6)]
+
+mns <- mns[3:4]
+sds <- sds[3:4]
+
+uplim <- mns+sds
+lowlim <- mns-sds
+
+plot(1,xlim=c(0.5,2.5),ylim=c(.95*min(lowlim),1.05*max(uplim)),xaxt='n',xlab="",
+ main='',
+ ylab=expression("Saplings per ha"),cex.lab=1.5
+)
+
+wd=0.4
+for(i in 1:2){
+ segments(structure[i],lowlim[i],structure[i],uplim[i])
+ segments(structure[i]+wd/2,lowlim[i],structure[i]-wd/2,lowlim[i])
+ segments(structure[i]+wd/2,uplim[i],structure[i]-wd/2,uplim[i])
+}
+
+points(structure,mns,cex=4,bg=cols,pch=21)
+text(structure,(uplim+4),c('A','A'))
+
+dev.off()