## Script to plot results from all model temporal runs
## Author: Andrew Tredennick
## Email: atredenn@gmail.com
## Date created: 10-09-2015
## Clear the workspace...
rm(list=ls())
####
#### Load libraries -----------------------------------------------------------
####
library(ggplot2)
library(plyr)
library(reshape2)
####
#### Read in model simulations, 1 by 1, and stack -----------------------------
####
setwd("../results/yearlyforecasts/")
all_files <- list.files()
torms <- grep("MEAN", all_files)
all_files <- all_files[-torms]
longdf <- data.frame(model=NA,
scenario=NA,
paramset=NA,
year=NA,
cover=NA)
for(i in 1:length(all_files)){
ttt <- readRDS(all_files[i])
ttt2 <- as.data.frame(apply(ttt, c(2,1), mean))
colnames(ttt2) <- paste0("params", 1:ncol(ttt2))
ttt2$year <- c(2011:2098)
ttt3 <- melt(ttt2, id.vars="year")
ttt3$model <- strsplit(all_files[i], "_")[[1]][1]
ttt3$scenario <- strsplit(all_files[i], "_")[[1]][2]
colnames(ttt3) <- c("year", "paramset", "cover", "model", "scenario")
tmpout <- ttt3[,c("model", "scenario", "paramset", "year", "cover")]
longdf <- rbind(longdf, tmpout)
print(paste(i, "out of", length(all_files)))
}
longdf <- longdf[2:nrow(longdf),]
longall <- longdf
longall$parammod <- paste(longall$param_set, longall$model)
# longdf <- subset(longall, model %in% c("bcc-csm1-1", "ccsm4"))
longdf <- longall
meandf <- ddply(longdf, .(year, scenario), summarise,
avgcover=mean(cover, na.rm=TRUE),
upcover=quantile(cover, 0.95, na.rm=TRUE),
locover=quantile(cover, 0.05, na.rm=TRUE))
####
#### Read in a format observation data ----------------------------------------
####
obs_data <- read.csv("../../data/wy_sagecover_subset_noNA.csv")
obs_data <- subset(obs_data, Year>1984) # subsets out NA CoverLag values
obs_agg <- ddply(obs_data, .(Year), summarise,
avgcover=mean(Cover),
upcover=quantile(Cover, 0.95),
locover=quantile(Cover, 0.05))
####
#### Make the plot ------------------------------------------------------------
####
g1 <- ggplot()+
geom_line(data=obs_agg, aes(x=Year, y=avgcover))+
# geom_ribbon(data=obs_agg,
# aes(x=Year, ymin=locover, ymax=upcover),
# alpha=0.5)+
geom_ribbon(data=meandf,
aes(x=year, ymin=locover, ymax=upcover, fill=scenario),
alpha=0.5)+
geom_line(data=meandf, aes(x=year, y=avgcover, color=scenario))+
geom_vline(aes(xintercept=2011))+
ylab("Mean sagebrush cover (%)")+
# guides(color=TRUE,fill=TRUE)+
scale_y_continuous(limits=c(0,30))+
scale_fill_manual(values=c("tan","coral","darkred"),
name="IPCC \nScenario",
labels=c("RCP 4.5", "RCP 6.0", "RCP 8.5"))+
scale_color_manual(values=c("tan","coral","darkred"),
name="IPCC \nScenario",
labels=c("RCP 4.5", "RCP 6.0", "RCP 8.5"))+
ggtitle("A) Long-term forecast")+
theme_bw()+
theme(legend.position = c(0.4, 0.8))
# ggsave("../temporal_forecast_wpois.png", width = 5, height = 4, dpi = 120)
####
#### PLOT A SHORT TERM FORECAST FROM ONE GCM
####
one_gcm <- subset(longdf, model == "miroc5" & year < 2022 & scenario=="rcp85")
g2 <- ggplot()+
geom_line(data=obs_agg, aes(x=Year, y=avgcover))+
geom_line(data=one_gcm, aes(x=year, y=cover, group=paramset), alpha=0.5, color="darkred")+
geom_vline(aes(xintercept=2011))+
ylab("Mean sagebrush cover (%)")+
guides(color=FALSE)+
scale_y_continuous(limits=c(0,30))+
ggtitle("B) Short-term forecast")+
theme_bw()
# ggsave("../short_term_miroc.png", width = 5, height = 4, dpi=120)
library(gridExtra)
tiff("../../figures/figure6.tiff", width = 5, height = 8, units = "in", res = 200)
gout <- grid.arrange(g1,g2)
dev.off()
## for presentations
library(cowplot)
gout <- plot_grid(g1,g2,ncol=2)
ggsave("../../results/temporal_forecasts_presentation.pdf",plot = gout, width = 8.5, height = 4)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.