# Coldwater Habitat Reconstruction Model - MGLP Data    
library(data.table)
library(mgcv)
library(graphics)
library(ggplot2)
library(gridExtra)
library(grid)
library(ggmap)
library(maps)
library(randomForest)
library(readr)
library(dplyr)
library(tidyr)
library(ggmosaic)
library(tidyverse)
library(magrittr)
library(mapproj)
library(reshape2)
library(cowplot)
#install.packages("googledrive")
#library(googledrive)
#set up color scheme
tier.colors=c("#0571b0", "#f4a582", "#ca0020")
myPal=c("#081d58","#1d91c0", "#c7e9b4", "#fed976", "#feb24c", "#ef6548", "#7f0000" )
#read in lat/long
mglp=read_csv("MGLP_Centroids_Lake.csv")

# Input data
# =======================================================

final.data=read.csv("mglp_lakes_mgmt_class_new.csv", header=TRUE)
final.data$class=factor(final.data$class, levels=c("Primary refuge","Primary protection", "Primary restoration", "Secondary refuge", "Secondary protection", "Secondary restoration", "Limited potential"))


#new tiers 

tier.1=13.6357
tier.2=23.17183

#how many lakes switch tiers because of climate change?

final.data=final.data %>% mutate(resilience.cc=case_when(
  tier=="Tier 1"& tier.cc=="Tier 1"~"Resilient Tier 1",
  tier=="Tier 1"& tier.cc!="Tier 1"~"At risk Tier 1",
  tier=="Tier 2"& tier.cc=="Tier 2"~"Resilient Tier 2",
  tier=="Tier 2"& tier.cc!="Tier 2"~"At risk Tier 2",
  tier=="Tier 3"& tier.cc=="Tier 3"~"Tier 3"
))
final.data$resilience.cc=factor(final.data$resilience.cc,levels=c("Resilient Tier 1", "At risk Tier 1","Resilient Tier 2","At risk Tier 2", "Tier 3") )




table(final.data$tier, final.data$resilience.cc)
table(final.data$tier)

table(final.data$tier.cc)


windows()
ggplot(data=final.data)+geom_segment(alpha=.5,aes(x=air_temp_july, xend=median.july, y=tdo3_est, yend=tdo3.cc_est, colour=resilience.cc), arrow = arrow(length = unit(0.1,"cm")))+theme_classic()+scale_colour_manual("Climate resilience", values=c("#253494", "#41b6c4", "#fee391", "#fe9929", "#fc4e2a"))+geom_hline(yintercept=tier.1, colour=tier.colors[1], lty=2, lwd=1)+geom_hline(yintercept=tier.2, colour=tier.colors[2], lty=2, lwd=1)+xlab(expression("Mean July air temperature("*degree*"C)"))+ylab(expression("Temperature at 3 mg/l of "*O[2]*" ("*degree*"C)"))
ggsave("tdo3_current_to_midcentury_resilience.png", height=6, width=6, units="in")

#######################


windows()
final.data%>%filter(tier!="Tier 3") %>% ggplot()+geom_point(aes(x=temp_change, y=climate.resilience, colour=class))+theme_bw()+geom_abline(intercept=0,slope=1)+scale_colour_manual("", values=myPal)+theme(panel.grid = element_blank(), legend.position = "bottom", strip.background = element_blank())+xlab(expression("July air temperature change ("*degree*"C)"))+ylab(expression("Climate resilience ("*degree*"C)"))+facet_wrap(~tier, ncol=1)
ggsave("climate_resilience_new_profiles_new.png", height=8, width=6.5, units="in")

#histogram

windows()
final.data%>%filter(tier!="Tier 3") %>%ggplot( aes(climate.resilience, fill=resilience.cc))+geom_histogram( colour="darkgrey",binwidth=0.5)+facet_wrap(~tier_current)+theme_bw()+xlab("Climate resilience")+ylab("Number of lakes")+scale_fill_manual("Climate resilience", values=c("#253494", "#41b6c4", "#fee391", "#fe9929", "#fc4e2a"))+xlab(expression("Climate resilience ("*degree*"C)"))+theme(panel.grid = element_blank(), legend.position = c(.2,.75), strip.background = element_blank())
ggsave("Climate_resilience_histogram.png", height=4, width =5, units="in")




#######################################
#for figure 2

##current vs. future with uncertainty
#plot a sample for easier viewing
windows()
plot1=ggplot(sample_n(final.data, size = 1500, replace = F, seed=4), aes(tdo3_est, tdo3.cc_est, colour=tier.cc))+geom_abline(slope=1, intercept=0)+geom_hline(yintercept=tier.1, colour=tier.colors[1])+geom_vline(xintercept=tier.1, colour=tier.colors[1])+geom_hline(yintercept=tier.2, colour=tier.colors[2])+geom_vline(xintercept=tier.2, colour=tier.colors[2])+geom_errorbar(aes(ymin=TDO3.cc.lower.PI,ymax=TDO3.cc.upper.PI), alpha=.05)+geom_errorbarh(aes(xmin=TDO3.current.lower.PI, xmax=TDO3.current.upper.PI), alpha=.05)+geom_point(alpha=.3)+  theme_classic()+xlab(expression(atop("Contemporary temperature", "at 3 mg/L of " *O[2]*" ("*degree*"C)")))+ylab(expression(atop("Mid-century temperature", "at 3 mg/L of " *O[2]*" ("*degree*"C)")))+scale_colour_manual("Tier", values=tier.colors)+theme_bw()+theme(legend.position="none", panel.grid = element_blank(), strip.background = element_blank(), legend.background = element_blank(), axis.text = element_text(size=12), axis.title = element_text(size=12))
print(plot1)
#ggsave("current_future_tdo3_with_PIs_new.png", height=6, width=6, units="in")


#histogram showing tier counts
windows()
plot2=ggplot(final.data)+geom_bar(aes(tier, fill=tier), alpha=.4, width=.5, colour="grey40")+scale_fill_manual("", values=tier.colors)+geom_bar(aes(tier.cc, fill=tier.cc), width=.2, colour="black")+theme_bw()+theme(legend.position = "none", panel.grid = element_blank(), strip.background = element_blank())+xlab("")+ylab("Number of lakes")
#ggsave("stat_tier_counts_present_future_new.png", height=3, width=8, units="in")


#by state
windows()
ggplot(final.data)+geom_bar(aes(tier, fill=tier), alpha=.4, width=.5, colour="grey40")+scale_fill_manual("", values=tier.colors)+geom_bar(aes(tier.cc, fill=tier.cc), width=.2, colour="black")+theme_bw()+theme(legend.position = "none", panel.grid = element_blank(), strip.background = element_blank())+xlab("")+ylab("Number of lakes")+facet_wrap(~STATE, scales="free_y", ncol=4)
ggsave("stat_tier_counts_present_future_by_state_new.png", height=3, width=8, units="in")

#############################################################


#map tdo3
all_states <- map_data("state")
states <- subset(all_states, region %in% c("indiana" , "minnesota", "michigan", "wisconsin", "iowa", "illinois", "south dakota", "north dakota") )

windows()
p <- ggplot()+ggtitle("Contemporary")
p <- p + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )
p=p+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")
p=p+geom_point( data=final.data, aes(Lon, Lat,  colour=tdo3_est), size=1.5)
p=p+scale_colour_distiller(name=expression(atop("Temperature at 3 mg/l", "of " *O[2]*" ("*degree*"C)")), palette="RdYlBu")
p=p+theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position=c(.2,.25) , legend.key = element_blank(),legend.background = element_blank(),
  legend.title=element_text(size=14),legend.text=element_text(size=14),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(),
  plot.title = element_text(face = "bold", size = (15), hjust=.5, vjust=-1),
)
p
ggsave("MGLP_tdo3_map_new.png", height=5, width=6, units="in")


#UCI
windows()
p <- ggplot()+ggtitle("Contemporary - UCI")
p <- p + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )
p=p+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")
p=p+geom_point( data=final.data, aes(Lon, Lat,  colour=tdo3_uci), size=1.5)
p=p+scale_colour_distiller(name=expression(atop("Temperature at 3 mg/l", "of " *O[2]*" ("*degree*"C)")), palette="RdYlBu")
p=p+theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position=c(.2,.25) , legend.key = element_blank(),legend.background = element_blank(),
  legend.title=element_text(size=14),legend.text=element_text(size=14),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(),
  plot.title = element_text(face = "bold", size = (15), hjust=.5, vjust=-1),
)
p
ggsave("MGLP_tdo3_UCI_map_new.png", height=5, width=6, units="in")

#LCI
windows()
p <- ggplot()+ggtitle("Contemporary - LCI")
p <- p + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )
p=p+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")
p=p+geom_point( data=final.data, aes(Lon, Lat,  colour=tdo3_lci), size=1.5)
p=p+scale_colour_distiller(name=expression(atop("Temperature at 3 mg/l", "of " *O[2]*" ("*degree*"C)")), palette="RdYlBu")
p=p+theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position=c(.2,.25) , legend.key = element_blank(),legend.background = element_blank(),
  legend.title=element_text(size=14),legend.text=element_text(size=14),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(),
  plot.title = element_text(face = "bold", size = (15), hjust=.5, vjust=-1),
)
p
ggsave("MGLP_tdo3_LCI_map_new.png", height=5, width=6, units="in")


#future TDO3

windows()
p <- ggplot()+ggtitle("Mid-century")
p <- p + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )
p=p+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")
p=p+geom_point( data=final.data, aes(Lon, Lat,  colour=tdo3.cc_est), size=1.5)
p=p+scale_colour_distiller(name=expression(atop("Temperature at 3 mg/l", "of " *O[2]*" ("*degree*"C)")), palette="RdYlBu")
p=p+theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position=c(.2,.25) , legend.key = element_blank(),legend.background = element_blank(),
  legend.title=element_text(size=14),legend.text=element_text(size=14),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(),
  plot.title = element_text(face = "bold", size = (15), hjust=.5, vjust=-1),
)
p
ggsave("MGLP_tdo3_future_map_new.png", height=5, width=6, units="in")

#map of Tiers with uncertainty
windows()
plot3 <- ggplot()+ggtitle("Contemporary") + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")+geom_point( data=final.data, aes(Lon, Lat,  colour=tier, alpha=tier.prob), size=1.5) +scale_colour_manual("", values=tier.colors)+ guides(alpha = FALSE) +theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position=c(.2,.25) , legend.key = element_blank(),legend.background = element_blank(),
  legend.title=element_blank(),legend.text=element_text(size=11),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(), 
  plot.title = element_text(face = "bold", size = (12), hjust=.5, vjust=-1),
)
plot3
#ggsave("MGLP_tier_map_contemporary_uncertainty_new.png", height=5, width=6, units="in")


#future climate
windows()
plot4 <- ggplot()+ggtitle("Mid-century")+ geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")+geom_point( data=final.data, aes(Lon, Lat,  colour=tier.cc, alpha=tier.prob.cc), size=1.5)+scale_colour_manual("Tier", values=tier.colors)+ theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position="none" , legend.key = element_blank(),
  #legend.title=element_text(size=16, hjust=-0.3),legend.text=element_text(size=16),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(),
  plot.title = element_text(face = "bold", size = (12), hjust=.5, vjust=-1),
)
plot4
#ggsave("MGLP_tier_map_midcentury_map_uncertainty_new.png", height=5, width=6, units="in")


##########################################################
#cowplot for multipanel figure
windows()
plot_grid(plot3, plot1,plot4,plot2, labels=c('A', 'B', 'C', 'D'), ncol=2,  rel_widths = c(1.3, 1))
ggsave("multipanel_tier_current_future_new.png", height=7, width=8, units = "in")
############################################
#calculate change in probability of being tier 1
final.data$prob.1.change=final.data$tier1.prob.cc-final.data$tier1.prob

#plot showing probability of being tier 1 as a function of both temp and watershed disturbance
#windows()
temp.prob.plot=ggplot(final.data, aes(median.july, tier1.prob.cc, colour=class))+geom_point(alpha=.7)+scale_colour_manual("",values=myPal)+theme_bw()+ylab("")+xlab(expression("Mid-century July air temperature ("*degree*"C)"))+theme(panel.grid = element_blank())

#windows()
ws.prob.plot=ggplot(final.data, aes(watershed_dist, tier1.prob.cc, colour=class))+geom_point(alpha=.7)+scale_colour_manual("",values=myPal)+theme_bw()+ylab("Probability Tier 1")+xlab("Proportion watershed disturbance")+theme(panel.grid = element_blank())

legend <- get_legend(
  # create some space to the left of the legend
  ws.prob.plot  +guides(color = guide_legend(nrow = 2)) +
    theme(legend.position = "top"))

##########################################################
#cowplot for multipanel figure
windows()
plots=plot_grid(ws.prob.plot+theme(legend.position = "none"), temp.prob.plot+theme(legend.position = "none"), labels=c('A', 'B'), ncol=2,  align='vh')
plot_grid(plots, legend, rel_heights =  c(1, .1), ncol=1)
ggsave("prob_tier_1_ws_temp.png", height=5, width=8, units = "in")
############################################





#map of probability of being Tier 1
#map of Tiers with uncertainty
windows()
tier1.prob.plot <- ggplot()+ggtitle("Contemporary") + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")+geom_point( data=filter(final.data, tier1.prob>.1), aes(Lon, Lat,   alpha=tier1.prob), size=1.5, colour=tier.colors[1]) +scale_alpha_continuous("Probability Tier 1")+theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position="none" , legend.key = element_blank(),legend.background = element_blank(),
  legend.title=element_text(size=11),legend.text=element_text(size=11),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(), 
  plot.title = element_text(face = "bold", size = (12), hjust=.5, vjust=-1),
)
tier1.prob.plot
#ggsave("MGLP_tier1_probability_map_contemporary_uncertainty_new.png", height=5, width=6, units="in")


#future climate
windows()
tier1.cc.prob.plot <- ggplot()+ggtitle("Future") + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")+geom_point( data=filter(final.data, tier1.prob.cc>.1), aes(Lon, Lat,   alpha=tier1.prob.cc), size=1.5, colour=tier.colors[1]) +scale_alpha_continuous("Probability Tier 1")+theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position="bottom" , legend.key = element_blank(),legend.background = element_blank(),
  legend.title=element_text(size=11),legend.text=element_text(size=11),
  axis.title.y=element_blank(),
  axis.title.x=element_blank(), 
  plot.title = element_text(face = "bold", size = (12), hjust=.5, vjust=-1),
)
tier1.cc.prob.plot
#ggsave("MGLP_tier_map_midcentury_map_uncertainty_new.png", height=5, width=6, units="in")

#histograms of probability
final.data$tier.1.prob.change=final.data$tier1.prob.cc-final.data$tier1.prob

windows()
current.tier1.hist=ggplot(data=filter(final.data, tier1.prob>.1),aes(tier1.prob))+geom_histogram(binwidth=.05, fill=tier.colors[1], colour="white")+theme_bw()+xlab("Tier 1 probability")+ylab("Lakes")+#ylim(lim=c(0,500))+
  theme(panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.title=element_text(size=11), 
  )
windows()
cc.tier1.hist=ggplot(data=filter(final.data, tier1.prob.cc>.1),aes(tier1.prob.cc))+geom_histogram(binwidth=.05, colour=tier.colors[1], fill=tier.colors[1], alpha=.2)+theme_bw()+xlab("Tier 1 probability")+ylab("Lakes")+ylim(lim=c(0,700))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        axis.title=element_text(size=11), 
  )


plot.with.inset=ggdraw()+draw_plot(tier1.prob.plot)+draw_plot(current.tier1.hist,x=-0.05, y=.08, width=.4, height=.35 )
plot.with.inset2=ggdraw()+draw_plot(tier1.cc.prob.plot)+draw_plot(cc.tier1.hist,x=-.05, y=.15, width=.4, height=.3 )

##########################################################
#cowplot for multipanel figure
windows()
plot_grid(plot.with.inset, plot.with.inset2,  labels=c('A', 'B'), ncol=1,  rel_widths = c(1, 1), rel_heights = c(1,1.05))
ggsave("multipanel_tier1_probs_current_future_new.png", height=6.5, width=4.5, units = "in")
############################################



####################################
#################################
##################################
#watershed resilience -

toplot=filter(final.data, !is.na(disturbance.change))

not.sensitive=anti_join(final.data, toplot)

#######################################
############



####combined figure with resilience and threshold

res.50=ggplot(toplot, aes(disturbance.change, fill=class))+geom_histogram(colour="black", binwidth=0.05)+facet_wrap(~tier.cc)+theme_bw()+theme(panel.grid = element_blank(), strip.background = element_blank(),strip.text = element_blank(), legend.position = "bottom", legend.background = element_blank())+scale_fill_manual("", values=myPal[c(2,3,5,6)])+xlab("Watershed resilience")+ylab("Lakes")

threshold.50=ggplot(toplot, aes(1-tier.threshold, fill=class))+geom_histogram(colour="black", binwidth=0.05)+facet_wrap(~tier.cc)+theme_bw()+theme(panel.grid = element_blank(), strip.background = element_blank(), legend.position = "none", legend.background = element_blank())+scale_fill_manual("", values=myPal[c(2,3,5,6)])+xlab("Watershed protection threshold")+ylab("Lakes")


#cowplot for multipanel figure
windows()
plot_grid(threshold.50, res.50, labels=c('A', 'B'), ncol=1)
ggsave("multipanel_watershed_protection_resilience50.png", height=6, width=8, units = "in")
############################################


###################################################

#map management tiers
#final.data=merge(final.data, mglp, by.x="mglp_id", by.y="LAKE_ID")


all_states <- map_data("state")
states <- subset(all_states, region %in% c("indiana" , "minnesota", "michigan", "wisconsin", "iowa", "illinois", "south dakota", "north dakota") )



windows()
p <- ggplot()
p <- p + geom_polygon( data=states, aes(x=long, y=lat, group = group),colour="black", fill="white" )
#p=p+geom_point( data=mglp, aes(Lon, Lat), size=.5, colour="grey90")
p=p+geom_point( data=final.data, aes(Lon, Lat,  colour=class, size=class))
#add back blue lakes for emphasis
p=p+geom_point( data=subset(final.data, class=="Primary protection"|class=="Primary refuge"|class=="Primary restoration"), aes(Lon, Lat,  colour=class, size=class))
p=p+scale_colour_manual("Management class",values =myPal )+scale_size_manual("Management class", values=c(2, 2, 2, 1.5,1.5, 1.5, 0.5))
p=p+theme(
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  axis.text= element_blank(),
  legend.position=c(.15,.3) , legend.key = element_blank(),
  #legend.title=element_text(size=16, hjust=-0.3),legend.text=element_text(size=16),
  axis.title.y=element_blank(),
  axis.title.x=element_blank()
)
p
ggsave("MGLP_new_management_class_map_all_lakes_new_profiles_new.png", height=6, width=8, units="in")




#histogram showing all lakes management class by state
windows()
ggplot(final.data)+geom_bar(aes(class, fill=class),  width=.5, colour="black")+scale_fill_manual("", values=myPal)+theme_bw()+theme(legend.position = "none", panel.grid = element_blank(), strip.background = element_blank(), axis.text.x=element_text(angle=90))+xlab("")+ylab("Number of lakes")+facet_wrap(~STATE, scales="free_y", ncol=2)
ggsave("management_class_counts_by_state_new.png", height=6, width=4, units="in")



##############################################################
#individual lake examples

#for individual lake figures, best to envision change in TDO3 per change in watershed development with uncertainty
scenario.plot.data=read.csv("ws_gradient_tdo3_new_PIs.csv")
ws.resilience=read_csv("tier_probs_watershed_gradient_new.csv")
climate.resilience=read_csv( "tier_probs_climate_gradient_new.csv")


scenario.plot.data.climate=read_csv( "temp_gradient_tdo3_new_PIs.csv")


#################################################

#make one with both and fewer lakes
test.lakes2=c( "WI600002303", "MN69033000","MIclntn46951" ,"MN29003600","MN38033000","IN155638231"  )


prob.plot2=ws.resilience[ws.resilience$mglp_id%in%test.lakes2,]
prob.plot2$tier.cc=paste("Tier", prob.plot2$tier, sep=" ")

plot.points2=final.data[final.data$mglp_id%in%test.lakes2,]

classes=select(final.data, mglp_id, class)
prob.plot2=merge(prob.plot2, classes)
prob.plot2$class=factor(prob.plot2$class, levels=c("Primary refuge",
                                                 "Primary protection", 
                                                 "Primary restoration",
                                                 "Secondary refuge", 
                                                 "Secondary protection",
                                                 "Secondary restoration" ))





p2=ggplot(prob.plot2, aes(ws_dist*100, prob))+geom_line(lwd=.5, aes( colour=tier.cc), alpha=.5)+theme_classic()+scale_colour_manual("Tier", values=tier.colors)+theme(legend.position = "right")+guides(colour = guide_legend(override.aes = list(alpha=1)))+xlab("Watershed disturbance")+ylab("Probability")+geom_point(data=plot.points2, aes(watershed_dist*100,tier.prob.cc , colour=tier.cc), size=4)+facet_wrap(class~mglp_id, nrow=1,dir="v")+theme(panel.background = element_rect(colour="grey"), strip.background = element_blank(), axis.text=element_text(size=10), axis.title = element_text(size=12), legend.position = "bottom",plot.margin = ggplot2::margin(b=0))+xlab("Watershed disturbance (%)")



p2

#prob with climate change
prob.plot3=climate.resilience[climate.resilience$mglp_id%in%test.lakes2,]
prob.plot3$tier=paste("Tier", prob.plot3$tier, sep=" ")


classes=select(final.data, mglp_id, class)
prob.plot3=merge(prob.plot3, classes)
prob.plot3$class=factor(prob.plot3$class, levels=c("Primary refuge",
                                                   "Primary protection", 
                                                   "Primary restoration",
                                                   "Secondary refuge", 
                                                   "Secondary protection",
                                                   "Secondary restoration" ))

#
plot.points2$class=factor(plot.points2$class, levels=c("Primary refuge",
                                                     "Primary protection", 
                                                     "Primary restoration",
                                                     "Secondary refuge", 
                                                     "Secondary protection",
                                                     "Secondary restoration" ))





p3=ggplot(prob.plot3, aes(july.temp, prob))+geom_line(lwd=.5, aes(colour=tier), alpha=.5)+theme_classic()+scale_colour_manual("Tier", values=tier.colors)+theme(legend.position = "right")+guides(colour = guide_legend(override.aes = list(alpha=1)))+xlab("July temperature")+ylab("Probability")+geom_point(data=plot.points2, aes(air_temp_july,tier.prob, colour=tier ), shape=21, fill="white",size=4)+geom_point(data=plot.points2, aes(median.july,tier.prob.cc, colour=tier.cc ), size=4)+facet_wrap(class~mglp_id, nrow=1)+theme(panel.background = element_rect(colour="grey"), strip.background = element_blank(),strip.text = element_blank(), axis.text=element_text(size=10), axis.title = element_text(size=12), legend.position = "bottom", plot.margin = ggplot2::margin(t=0))+xlab(expression("July air temperature ("*degree*"C)"))


legend <- get_legend(
  # create some space to the left of the legend
  p3  +guides(color = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom"))

##########################################################
#cowplot for multipanel figure
windows()
plots.probs=plot_grid(p2+theme(legend.position = "none"), p3+theme(legend.position = "none"), labels=c('A', 'B'), ncol=1,  align='vh')
plot_grid(plots.probs, legend, rel_heights =  c(1, .1), ncol=1)
ggsave("example_lakes_tdo3_tier_prob.png", height=6, width=8, units = "in")



###############################################



##################################################
#faceted table of numbers


count.table=data.frame(table(final.data$tier, final.data$tier.cc))
colnames(count.table)=c("tier", "tier.cc", "Counts")

count.table2=summarise(group_by(final.data, tier, tier.cc, class), counts=length(class))
count.table2$x=1
count.table2$y=c(.85, .35, .99, .5, .8, .2, .98, .4, .5)
count.table2$text.colour=c("white", "black","black","black","black","black","black","white","white" )



#what about a faceted bar chart or tiles
windows()
ggplot(final.data)+geom_bar(aes(1, fill=class), position="fill")+facet_grid(tier.cc~tier)+theme_bw()+theme(plot.margin = unit(c(1, 3, 1, 1), "lines"),axis.ticks=element_blank(), legend.position="bottom", panel.grid = element_blank(), panel.border = element_rect(colour="grey"),strip.background=element_blank(), axis.text = element_blank(), plot.title = element_text(hjust=0.5))+labs(x = "", y = "", title = "Current tier")+geom_text(data=count.table2, aes(x=x, y=y, label=counts, colour=text.colour))+scale_fill_manual("",values=myPal)  +scale_colour_manual(values=c("black", "white"))+guides(colour=FALSE)

grid.text("Mid-century tier", x = unit(0.95, "npc"), y = unit(0.50, "npc"), rot=-90)   
# ggsave("Cisco lake tiers grid with restoration actions.png", height=6, width=6, units="in")