##R code for analysis and map and figure creation for 
##Walters, A.W, Mandeville, C.P., and Rahel, F.R. The interaction of exposure and warming tolerance determines fish species vulnerability to ##warming stream temperatures. Biology Letters.

###########################
### Maps for CC MS ########
### last edit August 2 2018####
### Caitlin Mandeville ###
##########################

#### SET UP ####
library(maps)
library(maptools)
library(rgdal)
library(scales)
library(RColorBrewer)
library(raster)
library(sp)
library(rgeos)
library(car)
library(grDevices)
library(devEMF)

#read in river/stream data (obtained from http://www.nws.noaa.gov/geodata/catalog/hydro/html/rivers.htm in 2014)
#more recent versions available at https://www.weather.gov/gis/Rivers

Rivers.all<- readOGR(".", "rs14fe02")
Rivers<- Rivers.all[(Rivers.all$HUC2==14 | Rivers.all$HUC2==10 | Rivers.all$HUC2==16 | Rivers.all$HUC2==17),]
Streams.all<- readOGR(".", "rv14fe02")
Streams<- Streams.all[Streams.all$HUC > 10000000 & Streams.all$HUC <18000000,]

#### PREPARE FINAL (1559-SITE) DATABASE FOR MAPPING ####

# Enter database of 1559 sites (database exported from ArcMap after combining data sources, removing all sites w/ no fish data and all sites w/ #temp -9999), available at Dryad Digital Repository: https://doi.org/10.5061/dryad.76p9017

sites<- read.csv("sites1559.csv", stringsAsFactors=FALSE)
 
# CALCULATING LOWEST WARMING TOLERANCE AT EACH SITE
# 1. Redefine presence absence database such that, for all species, absences are 0 and presences are coded as that species' MWAT
sites$BBT<- recode(sites$BBT, "0 = NA ; 1 : hi = 19.6")
sites$BHS<- recode(sites$BHS, "0 = NA ; 1 : hi = 25.0")
sites$BKT<- recode(sites$BKT, "0 = NA ; 1 : hi = 18.3")
sites$BNT<- recode(sites$BNT, "0 = NA ; 1 : hi = 19.3")
sites$CCF<- recode(sites$CCF, "0 = NA ; 1 : hi = 32.4")
sites$CRP<- recode(sites$CRP, "0 = NA ; 1 : hi = 29.8")
sites$FHM<- recode(sites$FHM, "0 = NA ; 1 : hi = 28.8")
sites$FMS<- recode(sites$FMS, "0 = NA ; 1 : hi = 26.2")
sites$LSC<- recode(sites$LSC, "0 = NA ; 1 : hi = 23.8")
sites$MTS<- recode(sites$MTS, "0 = NA ; 1 : hi = 22.0")
sites$RBT<- recode(sites$RBT, "0 = NA ; 1 : hi = 19.4")
sites$SMB<- recode(sites$SMB, "0 = NA ; 1 : hi = 28.9")
sites$UTC<- recode(sites$UTC, "0 = NA ; 1 : hi = 26.0")
sites$WHS<- recode(sites$WHS, "0 = NA ; 1 : hi = 27.7")
sites$CUT<- recode(sites$CUT, "0 = NA ; 1 : hi = 18.1")
sites$CSH<- recode(sites$CSH, "0 = NA ; 1 : hi = 24.2")
sites$GSF<- recode(sites$GSF, "0 = NA ; 1 : hi = 31.1")
sites$GZS<- recode(sites$GZS, "0 = NA ; 1 : hi = 30.3")
sites$JDT<- recode(sites$JDT, "0 = NA ; 1 : hi = 25.4")
sites$NRH<- recode(sites$NRH, "0 = NA ; 1 : hi = 26.0")
sites$PMN<- recode(sites$PMN, "0 = NA ; 1 : hi = 33.0")
sites$RKB<- recode(sites$RKB, "0 = NA ; 1 : hi = 30.2")
sites$STR<- recode(sites$STR, "0 = NA ; 1 : hi = 27.8")
sites$YEP<- recode(sites$YEP, "0 = NA ; 1 : hi = 25.0")

# 2. Find lowest non-zero number in each row to get lowest MWAT at that site

sites$lowest.MWAT<- apply(sites[,1:24], 1, FUN=min, na.rm=TRUE)

# 3. Use this equation (sites$lowest.MWAT - S32_2080D) to get the lowest warming tolerance at each site

sites$lowest.WTOL<- (sites$lowest.MWAT - sites$S32_2080D)

# Bin sites$lowest.WTOL by rounding down to nearest integer
range(sites$lowest.WTOL) #warming tolerance ranges from -4.27 to 13.90 (= 19 bins)
sites$lowest.WTOL.bin<- floor(sites$lowest.WTOL)
table(sites$lowest.WTOL.bin)

hist(sites$lowest.WTOL, breaks=18) # compare histograms to visually double check binning
hist(sites$lowest.WTOL.bin, breaks=18)


# Bin modeled current temps by rounding down to nearest integer
range(sites$S1_93_11) #temps range from 4.43 to 23.12 (= 20 bins)
sites$TempBin<- .bincode(sites$S1_93_11, breaks=c(1:24), right=TRUE, include.lowest=FALSE)
hist(sites$S1_93_11, breaks=24) # compare histograms to visually double check binning
hist(sites$TempBin, breaks=24)

#### MAP ONE ####

# Base map
emf(file = "WYMapCurrentModTemps2_openpts.emf", width = 12, height = 9)
pdf("WYMapCurrentModTemps2_openpts.pdf", width=12, height=9)
par(mar=c(1,1,1,0))

map('state', xlim=c(-111,-104), ylim=c(41,45), fill=T, col="gray99", border="gray99")

plot(Streams, add=T, col="grey39")
plot(Rivers.all, add=T, col="grey10", lwd=1.8)

dev.off()

# Color palette
color.new<- colorRampPalette(c("darkblue","slategray3","red2"))
color.new1<- color.new(40)
palette(color.new1)

plot(1:40, 1:40, col=1:40, pch=19, cex=3) # view color spectrum
## For map, will use every other color for 0:10 and 20:30, and every color for 10:20, in order to use more neutral colors on the map

# Add sites
points(sites[sites$TempBin==4,]$X, sites[sites$TempBin==4,]$Y, col="#00008B", pch=16, cex=1.5)
points(sites[sites$TempBin==5,]$X, sites[sites$TempBin==5,]$Y, col="#181C95", pch=16, cex=1.5)
points(sites[sites$TempBin==6,]$X, sites[sites$TempBin==6,]$Y, col="#30389F", pch=16, cex=1.5)
points(sites[sites$TempBin==7,]$X, sites[sites$TempBin==7,]$Y, col="#414AA6", pch=16, cex=1.5)
points(sites[sites$TempBin==8,]$X, sites[sites$TempBin==8,]$Y, col="#515DAC", pch=16, cex=1.5)
points(sites[sites$TempBin==9,]$X, sites[sites$TempBin==9,]$Y, col="#6170B3", pch=16, cex=1.5)
points(sites[sites$TempBin==10,]$X, sites[sites$TempBin==10,]$Y, col="#7282BA", pch=16, cex=1.5)
points(sites[sites$TempBin==11,]$X, sites[sites$TempBin==11,]$Y, col="#8295C1", pch=16, cex=1.5)
points(sites[sites$TempBin==12,]$X, sites[sites$TempBin==12,]$Y, col="#92A8C7", pch=16, cex=1.5)
points(sites[sites$TempBin==13,]$X, sites[sites$TempBin==13,]$Y, col="#9AB1CB", pch=16, cex=1.5)
points(sites[sites$TempBin==14,]$X, sites[sites$TempBin==14,]$Y, col="#A1B1C7", pch=16, cex=1.5)
points(sites[sites$TempBin==15,]$X, sites[sites$TempBin==15,]$Y, col="#A5A8BD", pch=16, cex=1.5)
points(sites[sites$TempBin==16,]$X, sites[sites$TempBin==16,]$Y, col="#AD95A8", pch=16, cex=1.5)
points(sites[sites$TempBin==17,]$X, sites[sites$TempBin==17,]$Y, col="#B58293", pch=16, cex=1.5)
points(sites[sites$TempBin==18,]$X, sites[sites$TempBin==18,]$Y, col="#BD707E", pch=16, cex=1.5)
points(sites[sites$TempBin==19,]$X, sites[sites$TempBin==19,]$Y, col="#C55D69", pch=16, cex=1.5)
points(sites[sites$TempBin==20,]$X, sites[sites$TempBin==20,]$Y, col="#CD4A54", pch=16, cex=1.5)
points(sites[sites$TempBin==21,]$X, sites[sites$TempBin==21,]$Y, col="#D5383F", pch=16, cex=1.5)
points(sites[sites$TempBin==22,]$X, sites[sites$TempBin==22,]$Y, col="#E11C1F", pch=16, cex=1.5)
points(sites[sites$TempBin==23,]$X, sites[sites$TempBin==23,]$Y, col="#EE0000", pch=16, cex=1.5)

dev.off()

# Legend for color gradient
legend_image.new <- as.raster(c("#00008B","#181C95","#30389F","#414AA6","#515DAC","#6170B3","#7282BA",
                                "#8295C1","#92A8C7","#9AB1CB","#A1B1C7","#A5A8BD","#AD95A8","#B58293",
                                "#BD707E","#C55D69","#CD4A54","#D5383F","#E11C1F","#EE0000"), ncol=1)
plot(c(0,3),c(0,1),type = 'n', axes = F,xlab = '', ylab = '', main = '')
text(x=1.5, y = seq(1,5,l=5), labels = seq(1,5,l=5))
rasterImage(legend_image.new, 0,0,1,1)

dev.off()


#### MAP TWO ####

# Base map
emf(file = "WYMapLowestWarmingTolerance2_openpts.emf", width = 12, height = 9)
pdf("WYMapLowestWarmingTolerace.pdf", width=12, height=9)
par(mar=c(1,1,1,0))

map('state', xlim=c(-111,-104), ylim=c(41,45), fill=T, col="gray99", border="gray99")

plot(Streams, add=T, col="grey39")
plot(Rivers.all, add=T, col="grey10", lwd=1.8)

dev.off()

# Color palette
#Set separate color spectrum for negative and positive WTOL values
color.neg.wtol<- colorRampPalette(c("red2","navajowhite3"))
color.pos.wtol<- colorRampPalette(c("slategray2","darkblue"))

color.neg.wtol(7) # use 1:6 (leave out the final neutral color) for the neg/0 WTOL values
color.pos.wtol(14) # use 2:14 for pos WTOL values (leave out first neutral color)

# Add sites
points(sites[sites$lowest.WTOL.bin==0,]$X, sites[sites$lowest.WTOL.bin==0,]$Y, col="#D29573", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==1,]$X, sites[sites$lowest.WTOL.bin==1,]$Y, col="#AAC2E6", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==2,]$X, sites[sites$lowest.WTOL.bin==2,]$Y, col="#9CB2DE", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==3,]$X, sites[sites$lowest.WTOL.bin==3,]$Y, col="#8EA2D7", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==4,]$X, sites[sites$lowest.WTOL.bin==4,]$Y, col="#8092CF", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==5,]$X, sites[sites$lowest.WTOL.bin==5,]$Y, col="#7181C7", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==6,]$X, sites[sites$lowest.WTOL.bin==6,]$Y, col="#6371C0", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==7,]$X, sites[sites$lowest.WTOL.bin==7,]$Y, col="#5561B8", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==8,]$X, sites[sites$lowest.WTOL.bin==8,]$Y, col="#4751B1", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==9,]$X, sites[sites$lowest.WTOL.bin==9,]$Y, col="#3840A9", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==10,]$X, sites[sites$lowest.WTOL.bin==10,]$Y, col="#2A30A1", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==11,]$X, sites[sites$lowest.WTOL.bin==11,]$Y, col="#1C209A", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==12,]$X, sites[sites$lowest.WTOL.bin==12,]$Y, col="#0E1092", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==13,]$X, sites[sites$lowest.WTOL.bin==13,]$Y, col="#00008B", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==-5,]$X, sites[sites$lowest.WTOL.bin==-5,]$Y, col="#EE0000", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==-4,]$X, sites[sites$lowest.WTOL.bin==-4,]$Y, col="#E81D17", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==-3,]$X, sites[sites$lowest.WTOL.bin==-3,]$Y, col="#E33B2E", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==-2,]$X, sites[sites$lowest.WTOL.bin==-2,]$Y, col="#DD5945", pch=16, cex=1.5)
points(sites[sites$lowest.WTOL.bin==-1,]$X, sites[sites$lowest.WTOL.bin==-1,]$Y, col="#D8775C", pch=16, cex=1.5)

dev.off()

# Legend of color gradient

legend_image2<- as.raster(c(color.neg.wtol(7),color.pos.wtol(14)))
plot(c(0,3),c(0,1),type = 'n', axes = F,xlab = '', ylab = '', main = '')
text(x=1.5, y = seq(1,5,l=5), labels = seq(1,5,l=5))
rasterImage(legend_image2, 0,0,1,1)

dev.off()

#####################################
### Code for Table 1 and Figure 2 ###
### last edit July 31 2018 ####
#####################################

x<-read.csv("sites1559.csv") ##available at Dryad Digital Repository: https://doi.org/10.5061/dryad.76p9017

dim(x)
summary(x)
head(x)
tail(x)
str(x)

#differences between temperature scenarios
x$diff2040<-x$S30_2040D-x$S1_93_11
x$diff2080<-x$S32_2080D-x$S1_93_11
x$diff2040n<-x$S29_2040-x$S1_93_11
x$diff2080n<-x$S31_2080-x$S1_93_11

###################calculating warming tolerances

####SPECIES WITH MORE THAN 10 SITES
BBT<-x[x$BBT>0,]
dim(BBT)
BHS<-x[x$BHS>0,]
dim(BHS)
BKT<-x[x$BKT>0,]
dim(BKT)
BNT<-x[x$BNT>0,]
dim(BNT)
CCF<-x[x$CCF>0,]
dim(CCF)
CRP<-x[x$CRP>0,]
dim(CRP)
CSH<-x[x$CSH>0,]
dim(CSH)
CUT<-x[x$CUT>0,]
dim(CUT)
FHM<-x[x$FHM>0,]
dim(FHM)
FMS<-x[x$FMS>0,]
dim(FMS)
GSF<-x[x$GSF>0,]
dim(GSF)
GZS<-x[x$GZS>0,]
dim(GZS)
JDT<-x[x$JDT>0,]
dim(JDT)
LSC<-x[x$LSC>0,]
dim(LSC)
MTS<-x[x$MTS>0,]
dim(MTS)
NRH<-x[x$NRH>0,]
dim(NRH)
PMN<-x[x$PMN>0,]
dim(PMN)
RBT<-x[x$RBT>0,]
dim(RBT)
RKB<-x[x$RKB>0,]
dim(RKB)
SMB<-x[x$SMB>0,]
dim(SMB)
STR<-x[x$STR>0,]
dim(STR)
UTC<-x[x$UTC>0,]
dim(UTC)
WHS<-x[x$WHS>0,]
dim(WHS)
YEP<-x[x$YEP>0,]
dim(YEP)

BBT$WT<-19.59-BBT$S1_93_11
BBT$WT2040<-19.59-BBT$S30_2040D
BBT$WT2080<-19.59-BBT$S32_2080D
boxplot(BBT$WT, BBT$WT2040, BBT$WT2080)
mean(BBT$WT)
mean(BBT$WT2040)
mean(BBT$WT2080)
sd(BBT$WT)
sd(BBT$WT2040)
sd(BBT$WT2080)

BHS$WT<-25-BHS$S1_93_11
BHS$WT2040<-25-BHS$S30_2040D
BHS$WT2080<-25-BHS$S32_2080D
boxplot(BHS$WT, BHS$WT2040, BHS$WT2080)
mean(BHS$WT)
mean(BHS$WT2040)
mean(BHS$WT2080)
sd(BHS$WT)
sd(BHS$WT2040)
sd(BHS$WT2080)

BKT$WT<-18.34-BKT$S1_93_11
BKT$WT2040<-18.34-BKT$S30_2040D
BKT$WT2080<-18.34-BKT$S32_2080D
boxplot(BKT$WT, BKT$WT2040, BKT$WT2080)
mean(BKT$WT)
mean(BKT$WT2040)
mean(BKT$WT2080)
sd(BKT$WT)
sd(BKT$WT2040)
sd(BKT$WT2080)

BNT$WT<-19.32-BNT$S1_93_11
BNT$WT2040<-19.32-BNT$S30_2040D
BNT$WT2080<-19.32-BNT$S32_2080D
boxplot(BNT$WT, BNT$WT2040, BNT$WT2080)
mean(BNT$WT)
mean(BNT$WT2040)
mean(BNT$WT2080)
sd(BNT$WT)
sd(BNT$WT2040)
sd(BNT$WT2080)

CCF$WT<-32.38-CCF$S1_93_11
CCF$WT2040<-32.38-CCF$S30_2040D
CCF$WT2080<-32.38-CCF$S32_2080D
boxplot(CCF$WT, CCF$WT2040, CCF$WT2080)
mean(CCF$WT)
mean(CCF$WT2040)
mean(CCF$WT2080)
sd(CCF$WT)
sd(CCF$WT2040)
sd(CCF$WT2080)

CRP$WT<-29.84-CRP$S1_93_11
CRP$WT2040<-29.84-CRP$S30_2040D
CRP$WT2080<-29.84-CRP$S32_2080D
boxplot(CRP$WT, CRP$WT2040, CRP$WT2080)
mean(CRP$WT)
mean(CRP$WT2040)
mean(CRP$WT2080)
sd(CRP$WT)
sd(CRP$WT2040)
sd(CRP$WT2080)

CSH$WT<-24.2-CSH$S1_93_11
CSH$WT2040<-24.2-CSH$S30_2040D
CSH$WT2080<-24.2-CSH$S32_2080D
boxplot(CSH$WT, CSH$WT2040, CSH$WT2080)
mean(CSH$WT)
mean(CSH$WT2040)
mean(CSH$WT2080)
sd(CSH$WT)
sd(CSH$WT2040)
sd(CSH$WT2080)

CUT$WT<-18.1-CUT$S1_93_11
CUT$WT2040<-18.1-CUT$S30_2040D
CUT$WT2080<-18.1-CUT$S32_2080D
boxplot(CUT$WT, CUT$WT2040, CUT$WT2080)
mean(CUT$WT)
mean(CUT$WT2040)
mean(CUT$WT2080)
sd(CUT$WT)
sd(CUT$WT2040)
sd(CUT$WT2080)

FHM<-FHM[FHM$S1_93_11>0,]
FHM$WT<-28.78-FHM$S1_93_11
FHM$WT2040<-28.78-FHM$S30_2040D
FHM$WT2080<-28.78-FHM$S32_2080D
boxplot(FHM$WT, FHM$WT2040, FHM$WT2080)
mean(FHM$WT)
mean(FHM$WT2040)
mean(FHM$WT2080)
sd(FHM$WT)
sd(FHM$WT2040)
sd(FHM$WT2080)

FMS$WT<-26.22-FMS$S1_93_11
FMS$WT2040<-26.22-FMS$S30_2040D
FMS$WT2080<-26.22-FMS$S32_2080D
boxplot(FMS$WT, FMS$WT2040, FMS$WT2080)
mean(FMS$WT)
mean(FMS$WT2040)
mean(FMS$WT2080)
sd(FMS$WT)
sd(FMS$WT2040)
sd(FMS$WT2080)

GSF$WT<-31.05-GSF$S1_93_11
GSF$WT2040<-31.05-GSF$S30_2040D
GSF$WT2080<-31.05-GSF$S32_2080D
boxplot(GSF$WT, GSF$WT2040, GSF$WT2080)
mean(GSF$WT)
mean(GSF$WT2040)
mean(GSF$WT2080)
sd(GSF$WT)
sd(GSF$WT2040)
sd(GSF$WT2080)

GZS$WT<-30.3-GZS$S1_93_11
GZS$WT2040<-30.3-GZS$S30_2040D
GZS$WT2080<-30.3-GZS$S32_2080D
boxplot(GZS$WT, GZS$WT2040, GZS$WT2080)
mean(GZS$WT)
mean(GZS$WT2040)
mean(GZS$WT2080)
sd(GZS$WT)
sd(GZS$WT2040)
sd(GZS$WT2080)

JDT$WT<-25.37-JDT$S1_93_11
JDT$WT2040<-25.37-JDT$S30_2040D
JDT$WT2080<-25.37-JDT$S32_2080D
boxplot(JDT$WT, JDT$WT2040, JDT$WT2080)
mean(JDT$WT)
mean(JDT$WT2040)
mean(JDT$WT2080)
sd(JDT$WT)
sd(JDT$WT2040)
sd(JDT$WT2080)

LSC$WT<-23.79-LSC$S1_93_11
LSC$WT2040<-23.79-LSC$S30_2040D
LSC$WT2080<-23.79-LSC$S32_2080D
boxplot(LSC$WT, LSC$WT2040, LSC$WT2080)
mean(LSC$WT)
mean(LSC$WT2040)
mean(LSC$WT2080)
sd(LSC$WT)
sd(LSC$WT2040)
sd(LSC$WT2080)

MTS$WT<-21.95-MTS$S1_93_11
MTS$WT2040<-21.95-MTS$S30_2040D
MTS$WT2080<-21.95-MTS$S32_2080D
boxplot(MTS$WT, MTS$WT2040, MTS$WT2080)
mean(MTS$WT)
mean(MTS$WT2040)
mean(MTS$WT2080)
sd(MTS$WT)
sd(MTS$WT2040)
sd(MTS$WT2080)

NRH$WT<-26-NRH$S1_93_11
NRH$WT2040<-26-NRH$S30_2040D
NRH$WT2080<-26-NRH$S32_2080D
boxplot(NRH$WT, NRH$WT2040, NRH$WT2080)
mean(NRH$WT)
mean(NRH$WT2040)
mean(NRH$WT2080)
sd(NRH$WT)
sd(NRH$WT2040)
sd(NRH$WT2080)

PMN$WT<-32.94-PMN$S1_93_11
PMN$WT2040<-32.94-PMN$S30_2040D
PMN$WT2080<-32.94-PMN$S32_2080D
boxplot(PMN$WT, PMN$WT2040, PMN$WT2080)
mean(PMN$WT)
mean(PMN$WT2040)
mean(PMN$WT2080)
sd(PMN$WT)
sd(PMN$WT2040)
sd(PMN$WT2080)

RBT$WT<-19.4-RBT$S1_93_11
RBT$WT2040<-19.4-RBT$S30_2040D
RBT$WT2080<-19.4-RBT$S32_2080D
boxplot(RBT$WT, RBT$WT2040, RBT$WT2080)
mean(RBT$WT)
mean(RBT$WT2040)
mean(RBT$WT2080)
sd(RBT$WT)
sd(RBT$WT2040)
sd(RBT$WT2080)

RKB$WT<-30.2-RKB$S1_93_11
RKB$WT2040<-30.2-RKB$S30_2040D
RKB$WT2080<-30.2-RKB$S32_2080D
boxplot(RKB$WT, RKB$WT2040, RKB$WT2080)
mean(RKB$WT)
mean(RKB$WT2040)
mean(RKB$WT2080)
sd(RKB$WT)
sd(RKB$WT2040)
sd(RKB$WT2080)

SMB$WT<-28.9-SMB$S1_93_11
SMB$WT2040<-28.9-SMB$S30_2040D
SMB$WT2080<-28.9-SMB$S32_2080D
boxplot(SMB$WT, SMB$WT2040, SMB$WT2080)
mean(SMB$WT)
mean(SMB$WT2040)
mean(SMB$WT2080)
sd(SMB$WT)
sd(SMB$WT2040)
sd(SMB$WT2080)

STR$WT<-27.82-STR$S1_93_11
STR$WT2040<-27.82-STR$S30_2040D
STR$WT2080<-27.82-STR$S32_2080D
boxplot(STR$WT, STR$WT2040, STR$WT2080)
mean(STR$WT)
mean(STR$WT2040)
mean(STR$WT2080)
sd(STR$WT)
sd(STR$WT2040)
sd(STR$WT2080)

UTC$WT<-26.01-UTC$S1_93_11
UTC$WT2040<-26.01-UTC$S30_2040D
UTC$WT2080<-26.01-UTC$S32_2080D
boxplot(UTC$WT, UTC$WT2040, UTC$WT2080)
mean(UTC$WT)
mean(UTC$WT2040)
mean(UTC$WT2080)
sd(UTC$WT)
sd(UTC$WT2040)
sd(UTC$WT2080)

WHS<-WHS[WHS$S1_93_11>0,]
WHS$WT<-27.69-WHS$S1_93_11
WHS$WT2040<-27.69-WHS$S30_2040D
WHS$WT2080<-27.69-WHS$S32_2080D
boxplot(WHS$WT, WHS$WT2040, WHS$WT2080)
mean(WHS$WT)
mean(WHS$WT2040)
mean(WHS$WT2080)
sd(WHS$WT)
sd(WHS$WT2040)
sd(WHS$WT2080)

YEP<-YEP[YEP$S1_93_11>0,]
YEP$WT<-25.01-YEP$S1_93_11
YEP$WT2040<-25.01-YEP$S30_2040D
YEP$WT2080<-25.01-YEP$S32_2080D
boxplot(YEP$WT, YEP$WT2040, YEP$WT2080)
mean(YEP$WT)
mean(YEP$WT2040)
mean(YEP$WT2080)
sd(YEP$WT)
sd(YEP$WT2040)
sd(YEP$WT2080)

###############
###Figure 1 ###
###############

layout(matrix(c(1,2,3),3,1,byrow = TRUE), widths=c(1,1,1), heights=c(2,1,2))
layout.show(2)
par(mar=c(1,4,0.5,0.5), oma=c(11,1,0,0), xpd=NA)

boxplot(CUT$S1_93_11, BKT$S1_93_11,BNT$S1_93_11,RBT$S1_93_11, BBT$S1_93_11, MTS$S1_93_11, 
        LSC$S1_93_11,CSH$S1_93_11,BHS$S1_93_11,YEP$S1_93_11, JDT$S1_93_11, NRH$S1_93_11, UTC$S1_93_11, FMS$S1_93_11,
        WHS$S1_93_11, STR$S1_93_11,FHM$S1_93_11, SMB$S1_93_11,GZS$S1_93_11, CRP$S1_93_11, 
        RKB$S1_93_11, GSF$S1_93_11, CCF$S1_93_11,PMN$S1_93_11, axes=n, ylim=c(4,33),col="gray", ylab=expression(paste("Current temperature (",degree,"C)")),cex.lab=1.5, cex.axis=1.3)
points(1, 18.1, pch=17, col="black", add=TRUE, cex=1.5)
points(2, 18.3, pch=17, col="black", add=TRUE, cex=1.5)
points(3, 19.3, pch=17, col="black", add=TRUE, cex=1.5)
points(4, 19.4, pch=17, col="black", add=TRUE, cex=1.5)
points(5, 19.6, pch=17, col="black", add=TRUE, cex=1.5)
points(6, 22.0, pch=17, col="black", add=TRUE, cex=1.5)
points(7, 23.8, pch=17, col="black", add=TRUE, cex=1.5)
points(8, 24.2, pch=17, col="black", add=TRUE, cex=1.5)
points(9, 25.0, pch=17, col="black", add=TRUE, cex=1.5)
points(10, 25.0, pch=17, col="black", add=TRUE, cex=1.5)
points(11, 25.4, pch=17, col="black", add=TRUE, cex=1.5)
points(12, 26, pch=17, col="black", add=TRUE, cex=1.5)
points(13, 26, pch=17, col="black", add=TRUE, cex=1.5)
points(14, 26.2, pch=17, col="black", add=TRUE, cex=1.5)
points(15, 27.7, pch=17, col="black", add=TRUE, cex=1.5)
points(16, 27.8, pch=17, col="black", add=TRUE, cex=1.5)
points(17, 28.8, pch=17, col="black", add=TRUE, cex=1.5)
points(18, 28.9, pch=17, col="black", add=TRUE, cex=1.5)
points(19, 29.8, pch=17, col="black", add=TRUE, cex=1.5)
points(20, 30.2, pch=17, col="black", add=TRUE, cex=1.5)
points(21, 30.3, pch=17, col="black", add=TRUE, cex=1.5)
points(22, 31.1, pch=17, col="black", add=TRUE, cex=1.5)
points(23, 32.4, pch=17, col="black", add=TRUE, cex=1.5)
points(24, 33, pch=17, col="black", add=TRUE, cex=1.5)

text(.15,33,"A")
#adding exposure increases for 2080

boxplot(CUT$diff2080, BKT$diff2080, BNT$diff2080, RBT$diff2080, BBT$diff2080, 
        MTS$diff2080, LSC$diff2080, CSH$diff2080, BHS$diff2080, YEP$diff2080, JDT$diff2080,
        NRH$diff2080, UTC$diff2080, FMS$diff2080, STR$diff2080, WHS$diff2080, FHM$diff2080, 
        SMB$diff2080, CRP$diff2080, RKB$diff2080, GZS$diff2080, GSF$diff2080, 
        CCF$diff2080,PMN$diff2080,axes=n, ylab=expression(paste("Exposure by 2080 (",degree,"C) ")),  ylim=c(1.0,2.5), cex.lab=1.5,cex.axis=1.3)
text(.15,2.4,"B")


boxplot(CUT$WT2080, BKT$WT2080, BNT$WT2080, RBT$WT2080, BBT$WT2080, 
        MTS$WT2080, LSC$WT2080, CSH$WT2080, BHS$WT2080, YEP$diff2080, JDT$WT2080,
        NRH$WT2080, UTC$WT2080, FMS$WT2080, STR$WT2080, WHS$WT2080, FHM$WT2080, 
        SMB$WT2080, CRP$WT2080, RKB$WT2080, GZS$WT2080, GSF$WT2080, 
        CCF$WT2080,PMN$WT2080, axes=n, ylab=expression(paste("Warming tolerance by 2080 (",degree,"C)")), main="", ylim=c(-5.5,16), cex.lab=1.5,cex.axis=1.3)
axis(1, at=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21, 22, 23, 24), las=2,
     labels=c("Cutthroat trout", "Brook trout", "Brown trout", "Rainbow trout", "Burbot", 
              "Mountain sucker", "Leatherside chub", "Common shiner", "Bluehead sucker","Yellow perch", "Johnny darter", 
              "Shorthead redhorse","Utah chub","Flannelmouth sucker","Central stoneroller","White sucker",
              "Fathead minnow","Smallmouth bass", "Common carp", "Rock bass","Gizzard shad", 
              "Green sunfish","Channel catfish", "Plains minnow"), cex.axis=1.3)
segments(0,0,24,0, lty=2)
text(1,-5,"15.6%", cex=.8)
text(2,-5,"2.6%", cex=.8)
text(3,-5,"33.1%", cex=.8)
text(4,-5,"32.2%", cex=.8)
text(5,-5,"66.7%", cex=.8)
text(6,-5,"5.5%", cex=.8)
text(.15,16,"C")

########################################################
##mean exposure

mean(BBT$diff2040)
min(BBT$diff2040)
max(BBT$diff2040)
mean(BBT$diff2080)
min(BBT$diff2080)
max(BBT$diff2080)

mean(BHS$diff2040)
mean(BHS$diff2080)

mean(BKT$diff2040)
mean(BKT$diff2080)

mean(BNT$diff2040)
mean(BNT$diff2080)

mean(CCF$diff2040)
mean(CCF$diff2080)

mean(CRP$diff2040)
mean(CRP$diff2080)

mean(CSH$diff2040)
mean(CSH$diff2080)

mean(CUT$diff2040)
min(CUT$diff2040)
max(CUT$diff2040)
median(CUT$diff2040)
mean(CUT$diff2080)
min(CUT$diff2080)
max(CUT$diff2080)

mean(FHM$diff2040)
mean(FHM$diff2080)

mean(FMS$diff2040)
mean(FMS$diff2080)

mean(GSF$diff2040)
mean(GSF$diff2080)

mean(GZS$diff2040)
mean(GZS$diff2080)

mean(JDT$diff2040)
mean(JDT$diff2080)

mean(LSC$diff2040)
mean(LSC$diff2080)

mean(MTS$diff2040)
mean(MTS$diff2080)

mean(NRH$diff2040)
mean(NRH$diff2080)

mean(PMN$diff2040)
min(PMN$diff2040)
max(PMN$diff2040)
mean(PMN$diff2080)
min(PMN$diff2080)
max(PMN$diff2080)

mean(RBT$diff2040)
mean(RBT$diff2080)

mean(RKB$diff2040)
mean(RKB$diff2080)

mean(SMB$diff2040)
mean(SMB$diff2080)

mean(STR$diff2040)
mean(STR$diff2080)

mean(UTC$diff2040)
mean(UTC$diff2080)

mean(WHS$diff2040)
mean(WHS$diff2080)

mean(YEP$diff2040)
mean(YEP$diff2080)