5_Supplementary_figures_code.Rmd

---
title: "Supplementary Figures"
author: "SS"
date: "February 02, 2019"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r load_libraries}
library(MASS)
library(audio)
library(sp)
library(foreign)
library(rgdal)
library(maptools)
library(rgeos)
library(doParallel)
library(rasterVis)
library(dismo)
library(plotKML)
library(SDMTools)
library(PBSmapping)
library(lme4)
library(blme)
library(raster)
library(fields)
library(RColorBrewer)
library(sjmisc)
library(ncdf4)
library(knitr)
library(lattice)
library(mgcv)
library(coefplot)
library(R2jags)
library(ggplot2)
library(stringr)
library(GISTools)
```

main_working_directory="C:/Users/Shannon/Desktop/Ecoregions"
home<-"C:/Users/Shannon/Desktop/Ecoregions"

graphs_directory="C:/Users/Shannon/Desktop/Ecoregions/output"

ECO <- readOGR(file.path(home,'shapefiles','ecoregion_exportPolygon.shp')) # ecoregions
wlrd.p <- readOGR(file.path(home,'shapefiles','TM_WORLD_BORDERS_SIMPL_PC150.shp'))


```{r get rid of the holes in the polygon shapefiles}
ecos_list<-c()
for (i in 1:150){
  eco_i<-Polygons((Filter(function(f){f@ringDir==1}, ECO@polygons[[i]]@Polygons)), ID=i)
  ecos_list<-append(ecos_list, values=eco_i, after = length(ecos_list))
  #include a brief pause because if running in Rstudio, it takes a while for the code to run and for the value to be loaded into the global environment. If there is no pause, the next iteration of the loop starts before the previous value is fully saved and loaded into the environment, and there can be errors in the shapefile 
  Sys.sleep(.2)
}
ecos<-SpatialPolygons(ecos_list)

ecos$ERG<-ECO$ERG
ecos$Ecoregion<-ECO$Ecoregion
ecos@proj4string<-ECO@proj4string
ecos@plotOrder<-ECO@plotOrder
ecos@data<-ECO@data

ECO<-ecos
ecor.sg<-ECO
```

setwd("C:/Users/Shannon/Desktop/Ecoregions")
source(file = "HighstatLibV10.R")  
source(file = "MCMCSupportHighstatV4.R")
source(file= "MyBUGSOutput.R")
setwd(main_working_directory)

#Reef_Check_Bleaching_Data is all the Reef Check data. We filter out a few data points that are missing lat or lon, but we do not filter out the data that falls outside of ecoregions. This is used for the lat and long bleaching severity histograms.
Reef_Check_Bleaching_Data <- read.csv(file="Reef Check Data Raw.csv", header=TRUE, sep=",")
coral_diversity<-read.csv("coral_diversity.csv", header=TRUE, sep=",")
names(coral_diversity)[1]<-"Ecoregion"
setwd("C:/Users/Shannon/Desktop/Ecoregions")
mean_variables_per_ecoregion<- read.csv("EcoRegions_mean_variables.csv")
mean_variables_per_ecoregion<-subset(mean_variables_per_ecoregion, samples>=10)

#Bleaching_Data_with_cortad_variables filters out anything that is outside of the ecoregions. this is used for most of the other plots
StudyTitle<-"Reef_Check"
csv_Title<-paste(StudyTitle, "_with_cortad_variables_with_annual_rate_of_SST_change.csv", sep="")
Bleaching_Data_with_cortad_variables<-read.csv(csv_Title, header=TRUE, sep=",")
names(Bleaching_Data_with_cortad_variables)[names(Bleaching_Data_with_cortad_variables)=="ï..Reef.ID"]<-"Reef.ID"

```{r plotmap_function}
plot.map<- function(database,center,transf=T,...){
  Obj <- map(database,...,plot=F)
  coord <- cbind(Obj[[1]],Obj[[2]])
  newproj <- "+proj=merc +lon_0=150 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" #utm
  nextproj<-"+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" #latlong
  # split up the coordinates
  id <- rle(!is.na(coord[,1]))
  id <- matrix(c(1,cumsum(id$lengths)),ncol=2,byrow=T)
  polygons <- apply(id,1,function(i){coord[i[1]:i[2],]})
  
  # split up polygons that differ too much
  polygons <- lapply(polygons,function(x){
    x[,1] <- x[,1] + center
    x[,1] <- ifelse(x[,1]>180,x[,1]-360,x[,1])
    if(sum(diff(x[,1])>300,na.rm=T) >0){
      id <- x[,1] < 0
      x <- rbind(x[id,],c(NA,NA),x[!id,])
    }
    x
  })
  # reconstruct the object
  polygons <- do.call(rbind,polygons)
  
  
  colnames(polygons)<-c("x",'y')
  polygons<-as.data.frame(polygons)
  z<-complete.cases(polygons)
  p<-z
  z<-cbind(z,z)
  polygons<-polygons[complete.cases(polygons),]
  coordinates(polygons)<-~x+y
  proj4string(polygons)<-CRS(nextproj)
  if(transf==T){ polygons<-spTransform(polygons,CRS(newproj))}
  
  z[p==F,]<-c(NA,NA)
  z[which(p==T),]<-coordinates(polygons)
  Obj[[1]] <- z[,1]
  Obj[[2]] <- z[,2]
  
  map(Obj,...)
}
```

#format the data
```{r format data}
#format longitude
for (i in 1:length(Reef_Check_Bleaching_Data$Longitude.Degrees)){
  cell<-as.double(Reef_Check_Bleaching_Data$Longitude.Degrees[i])+ as.double(Reef_Check_Bleaching_Data$Longitude.Minutes[i]/60)+ as.double(Reef_Check_Bleaching_Data$Longitude.Seconds[i]/60/60)
  if (Reef_Check_Bleaching_Data$Longitude.Cardinal.Direction[i]=='W'){cell<-cell*-1}
  Reef_Check_Bleaching_Data$Longitude.Degrees[i]<-cell
}

#format latitude
for (i in 1:length(Reef_Check_Bleaching_Data$Latitude.Degrees)){
  cell<-as.double(Reef_Check_Bleaching_Data$Latitude.Degrees[i])+ as.double(Reef_Check_Bleaching_Data$Latitude.Minutes[i]/60)+ as.double(Reef_Check_Bleaching_Data$Latitude.Seconds[i]/60/60)
  if (Reef_Check_Bleaching_Data$Latitude.Cardinal.Direction[i]=='S'){cell<-cell*-1}
  Reef_Check_Bleaching_Data$Latitude.Degrees[i]<-cell
}

Reef_Check_Bleaching_Data = subset(Reef_Check_Bleaching_Data, select = -c(Latitude.Minutes,Latitude.Seconds,Latitude.Cardinal.Direction, Longitude.Minutes, Longitude.Seconds, Longitude.Cardinal.Direction))

#calculate average bleaching
average_bleaching<-array(data=NA,dim=length(Reef_Check_Bleaching_Data$S1))
for (i in 1:length(average_bleaching)){
num.transects<-0
bleaching_sum<-0
if (is.na(Reef_Check_Bleaching_Data$S1[i])==FALSE){num.transects=num.transects+1; bleaching_sum<-bleaching_sum+Reef_Check_Bleaching_Data$S1[i]}
if (is.na(Reef_Check_Bleaching_Data$S2[i])==FALSE){num.transects=num.transects+1; bleaching_sum<-bleaching_sum+Reef_Check_Bleaching_Data$S2[i]}
if (is.na(Reef_Check_Bleaching_Data$S3[i])==FALSE){num.transects=num.transects+1; bleaching_sum<-bleaching_sum+Reef_Check_Bleaching_Data$S3[i]}
if (is.na(Reef_Check_Bleaching_Data$S4[i])==FALSE){num.transects=num.transects+1; bleaching_sum<-bleaching_sum+Reef_Check_Bleaching_Data$S4[i]}
average_bleaching[i]<-bleaching_sum/num.transects
}
Reef_Check_Bleaching_Data$Average_bleaching<-as.double(average_bleaching)

#use only the population bleaching data, not colony data. remove data pts that have no lat or lon
population_rows<-which(Reef_Check_Bleaching_Data$Organism.Code=="Bleaching (% of population)")
Reef_Check_Bleaching_Data<-data.frame(Reef_Check_Bleaching_Data[population_rows,])
Reef_Check_Bleaching_Data<-subset(Reef_Check_Bleaching_Data, !is.na(Longitude.Degrees))
Reef_Check_Bleaching_Data<-subset(Reef_Check_Bleaching_Data, !is.na(Latitude.Degrees))

#fit data into ecoregions
library(raster)
 coordinates(Reef_Check_Bleaching_Data)<- ~Longitude.Degrees+Latitude.Degrees
 proj4string(Reef_Check_Bleaching_Data)<-"+proj=longlat +ellps=WGS84 +datum=WGS84"
 bldata<-spTransform(Reef_Check_Bleaching_Data,proj4string(ECO))
test<-over(bldata,ECO)
Bleaching_Data_with_cortad_variables$Region<-as.character(test$ERG)

ECO$num_studies<-0
for (i in 1:150){
  erg_code<-ECO$ERG[i]
  ECO$num_studies[i]<-table(Bleaching_Data_with_cortad_variables$Region)[as.character(erg_code)]
}

standardize_function<-function(x){
  x.standardized=(x-mean(na.omit(x)))/sd(na.omit(x))
  return(x.standardized)
}

#variables_to_standardize<-c("Temperature_Kelvin", "Depth", "Latitude.Degrees", "Longitude.Degrees", "Windspeed", "SSTA", "SSTA_Maximum", "SSTA_DHW", "SSTA_DHWMax", "SSTA_DHWMean", "SSTA_Frequency", "SSTA_FrequencyMax", "SSTA_FrequencyMean", "TSA", "TSA_Maximum", "TSA_Mean", "TsA_Frequency", "TSA_FrequencyMax", "TSA_FrequencyMean", "TSA_DHW", "TSA_DHWMax", "TSA_DHWMean")

#Latitude
#Bleaching_Data_with_cortad_variables$Latitude<-Bleaching_Data_with_cortad_variables$Latitude.Degrees
Bleaching_Data_with_cortad_variables$Latitude<-abs(Bleaching_Data_with_cortad_variables$Latitude.Degrees)
Bleaching_Data_with_cortad_variables$Latitude.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$Latitude)

#ClimSST
Bleaching_Data_with_cortad_variables$ClimSST.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$ClimSST)

#Temperature
Bleaching_Data_with_cortad_variables$Temperature_Kelvin.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$Temperature_Kelvin)

#Temperature_Mean
Bleaching_Data_with_cortad_variables$Temperature_Mean.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$Temperature_Mean)

#Temperature_Minimum
Bleaching_Data_with_cortad_variables$Temperature_Minimum.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$Temperature_Minimum)

#Temperature_Maximum
Bleaching_Data_with_cortad_variables$Temperature_Maximum.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$Temperature_Maximum)

#Depth
Bleaching_Data_with_cortad_variables$Depth.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$Depth)

#Year
Bleaching_Data_with_cortad_variables$Year.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$Year)

#Latitude
Bleaching_Data_with_cortad_variables$Latitude.Degrees.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$Latitude.Degrees)

#Longitude
Bleaching_Data_with_cortad_variables$Longitude.Degrees.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$Longitude.Degrees)

#WindSpeed
Bleaching_Data_with_cortad_variables$Windspeed.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$Windspeed)

#SSTA
Bleaching_Data_with_cortad_variables$SSTA.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA)

#SSTA_Minimum
Bleaching_Data_with_cortad_variables$SSTA_Minimum.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Minimum)

#SSTA_Maximum
Bleaching_Data_with_cortad_variables$SSTA_Maximum.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Maximum)

#SSTA_DHW
Bleaching_Data_with_cortad_variables$SSTA_DHW.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$SSTA_DHW)

#SSTA_DHWMax
Bleaching_Data_with_cortad_variables$SSTA_DHWMax.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$SSTA_DHWMax)

#SSTA_DHWMean
Bleaching_Data_with_cortad_variables$SSTA_DHWMean.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$SSTA_DHWMean)

#SSTA_Frequency
Bleaching_Data_with_cortad_variables$SSTA_Frequency.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Frequency)

#SSTA_FrequencyMax
Bleaching_Data_with_cortad_variables$SSTA_FrequencyMax.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_FrequencyMax)

#SSTA_FrequencyMean
Bleaching_Data_with_cortad_variables$SSTA_FrequencyMean.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_FrequencyMean)

#TSA
Bleaching_Data_with_cortad_variables$TSA.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA)

#TSA_Minimum
Bleaching_Data_with_cortad_variables$TSA_Minimum.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_Minimum)

#TSA_Maximum
Bleaching_Data_with_cortad_variables$TSA_Maximum.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_Maximum)

#TSA_Mean
Bleaching_Data_with_cortad_variables$TSA_Mean.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_Mean)

#TSA_Frequency
Bleaching_Data_with_cortad_variables$TSA_Frequency.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_Frequency)

#TSA_FrequencyMax
Bleaching_Data_with_cortad_variables$TSA_FrequencyMax.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_FrequencyMax)

#TSA_FrequencyMean
Bleaching_Data_with_cortad_variables$TSA_FrequencyMean.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_FrequencyMean)

#TSA_DHW
Bleaching_Data_with_cortad_variables$TSA_DHW.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_DHW)

#TSA_DHWMax
Bleaching_Data_with_cortad_variables$TSA_DHWMax.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_DHWMax)

#TSA_DHWMean
Bleaching_Data_with_cortad_variables$TSA_DHWMean.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$TSA_DHWMean)

#Temperature_Kelvin_Standard_Deviation
Bleaching_Data_with_cortad_variables$Temperature_Kelvin_Standard_Deviation.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$Temperature_Kelvin_Standard_Deviation)

#SSTA_Standard_Deviation
Bleaching_Data_with_cortad_variables$SSTA_Standard_Deviation.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Standard_Deviation)

#SSTA_DHW_Standard_Deviation
Bleaching_Data_with_cortad_variables$SSTA_DHW_Standard_Deviation.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Standard_Deviation)

#SSTA_Frequency_Standard_Deviation
Bleaching_Data_with_cortad_variables$SSTA_Frequency_Standard_Deviation.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Standard_Deviation)

#TSA_Standard_Deviation
Bleaching_Data_with_cortad_variables$TSA_Standard_Deviation.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Standard_Deviation)

#TSA_Frequency_Standard_Deviation
Bleaching_Data_with_cortad_variables$TSA_Frequency_Standard_Deviation.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Standard_Deviation)

#TSA_DHW_Standard_Deviation
Bleaching_Data_with_cortad_variables$TSA_DHW_Standard_Deviation.standardized<- standardize_function(Bleaching_Data_with_cortad_variables$SSTA_Standard_Deviation)

#rate_of_SST_change
Bleaching_Data_with_cortad_variables$rate_of_SST_change.standardized<-standardize_function(Bleaching_Data_with_cortad_variables$rate_of_SST_change)


Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(Temperature_Kelvin))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(Windspeed))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(SSTA))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(SSTA_Mean))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(SSTA_Maximum))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(SSTA_DHW))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(SSTA_DHWMax))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(Average_bleaching))
Bleaching_Data_with_cortad_variables <- subset(Bleaching_Data_with_cortad_variables, !is.na(Longitude.Degrees))

```

#Fig S1
```{r plot number of surveys per ecoregion}
setwd(graphs_directory)
D=ECO$num_studies
brks<-c(0,1,100,200,300,400,500,600,700,1000)
cols <- c("white", brewer.pal(9, "Blues")[2],brewer.pal(9, "Blues")[3:9])
grps<-as.ordered(cut(D, brks, include.lowest=TRUE))

tiff(file='Figure_S1_Final.tif',height=800,width=3000,res=300)
par(mgp=c(0.5,0.6,0), mar=c(1,1,1,1))
plot.map("world", center=0 , ylim=c(-5500000,5500000), fill=T, col="grey90",border='grey70',mar=c(4,4,1,18)) #center is still 0
axis(1,at=c(-10018754.17,3339584.724,16697920),lab=c('60°','180°','-60°'),las=1,tcl=0.35,mgp=c(-1,-1.3,0))
axis(2, at=c(23*111319.4666666667,0,-23*111319.4666666667),labels=c('23°','0°','-23°'),las=3,tcl=0.35,mgp=c(-2,-1.3,0),hadj=.4)
axis(3,at=c(-10018754.17,3339584.724,16697920),lab=c('','',''),las=1,tcl=0.35,mgp=c(-1,-1.3,0))
axis(4, at=c(23*111319.4666666667,0,-23*111319.4666666667),labels=c('','',''),las=2,tcl=0.35,mgp=c(-1,-0.6,0),hadj=0)
box()
plot(ECO, ylim=c(-5500000,5500000), col=cols[unclass(grps)],add=T,lwd=.5)
windowsFonts(A=windowsFont('Arial Unicode MS'))
text(-7868896,-2922012,'Indian Ocean',cex=0.7,family='A')
text(9438742,487176,'Pacific Ocean',cex=0.7,family='A')
plotrix::color.legend(9684797.171+25e5,-28*111319.4666666667,15807371.62+25e5,-23.5*111319.4666666667,legend=c(round(min(brks), digits=1),round(max(brks), digits=1)),rect.col=cols,cex=1)
north.arrow(x=(-16654136+111319.4*320), y=1615153*2, len=(111319.4*2), lab="N", cex=.7)
dev.off()

setwd(main_working_directory)
```

#Fig S2
```{r sites at lat and lon}
#number of sites at latitude
sites_at_lat<-cbind(Bleaching_Data_with_cortad_variables$Reef.ID, Bleaching_Data_with_cortad_variables$Latitude.Degrees)
for (i in 1:(dim(sites_at_lat)[1])){
sites_at_lat[i]<-paste(sites_at_lat[i,1], "LAT", sites_at_lat[i,2], sep = "")
}
sites_at_lat<-sites_at_lat[,1]
sites_at_lat<-unique(sites_at_lat)
for (i in 1:length(sites_at_lat)){
sites_at_lat[i]<-str_split(sites_at_lat[i], "LAT")[[1]][2]
}
sites_at_lat<-as.numeric(sites_at_lat)

setwd(graphs_directory)
steps_lat<-5
break_pts_lat<-seq(floor(min(na.omit(Bleaching_Data_with_cortad_variables$Latitude.Degrees)))-.5, ceiling(max(na.omit(Bleaching_Data_with_cortad_variables$Latitude.Degrees)))+steps_lat+.5, steps_lat)
tiff("Figure_S2_Final_Latitude.tif",res=300,width=1600,height=1200)
hist(sites_at_lat, xlab="Latitude", ylab="Number of Sites", breaks=break_pts_lat, main="")
dev.off()
setwd(main_working_directory)

#number of sites at longitude
sites_at_lon<-cbind(Bleaching_Data_with_cortad_variables$Reef.ID, Bleaching_Data_with_cortad_variables$Longitude.Degrees)
for (i in 1:(dim(sites_at_lon)[1])){
sites_at_lon[i]<-sites_at_lon[i]<-paste(sites_at_lon[i,1], "LON", sites_at_lon[i,2], sep = "")
}
sites_at_lon<-sites_at_lon[,1]
sites_at_lon<-unique(sites_at_lon)
for (i in 1:length(sites_at_lon)){
sites_at_lon[i]<-str_split(sites_at_lon[i], "LON")[[1]][2]
}
sites_at_lon<-as.numeric(sites_at_lon)

setwd(graphs_directory)
steps_lon<-5
break_pts_lon<-seq(floor(min(na.omit(Bleaching_Data_with_cortad_variables$Longitude.Degrees)))-.5, ceiling(max(na.omit(Bleaching_Data_with_cortad_variables$Longitude.Degrees)))+steps_lon+.5, steps_lon)
tiff("Figure_S2_Final_Longitude.tif",res=300,width=1600,height=1200)
hist(sites_at_lon, xlab="Longitude", ylab="Number of Sites", breaks=break_pts_lon, main="")
dev.off()
setwd(main_working_directory)
```

#Fig S3-S4
```{r plot severity bleaching at lat and lon}
Severity<-array(data=NA, dim=length(Bleaching_Data_with_cortad_variables$Average_bleaching))
for (i in 1:length(Bleaching_Data_with_cortad_variables$Average_bleaching)){
if(Bleaching_Data_with_cortad_variables$Average_bleaching[i]<1){Severity[i]=0}
if((Bleaching_Data_with_cortad_variables$Average_bleaching[i]>=1) & (Bleaching_Data_with_cortad_variables$Average_bleaching[i]<10)){Severity[i]=1}
if((Bleaching_Data_with_cortad_variables$Average_bleaching[i]>=10) & (Bleaching_Data_with_cortad_variables$Average_bleaching[i]<50)){Severity[i]=2}
if(Bleaching_Data_with_cortad_variables$Average_bleaching[i]>=50){Severity[i]=3}
}
Bleaching_Data_with_cortad_variables$Severity<-Severity

Bleaching_Data_Reef_Check_all= subset(Bleaching_Data_with_cortad_variables)
Bleaching_Data_Reef_Check_None=subset(Bleaching_Data_with_cortad_variables)
Bleaching_Data_Reef_Check_Mild<-subset(Bleaching_Data_with_cortad_variables, Severity>0)
Bleaching_Data_Reef_Check_Moderate<-subset(Bleaching_Data_with_cortad_variables, Severity>(1))
Bleaching_Data_Reef_Check_Severe<-subset(Bleaching_Data_with_cortad_variables, Severity>(2))
number_of_studies_Reef_Check=length(Bleaching_Data_with_cortad_variables$Severity)

break_pts_lat<-seq(-40, 40, 5)
lat0<-hist(na.omit(Bleaching_Data_Reef_Check_None$Latitude.Degrees), breaks=break_pts_lat)
lat1<-hist(na.omit(Bleaching_Data_Reef_Check_Mild$Latitude.Degrees), breaks=break_pts_lat)
lat2<-hist(na.omit(Bleaching_Data_Reef_Check_Moderate$Latitude.Degrees), breaks=break_pts_lat)
lat3<-hist(na.omit(Bleaching_Data_Reef_Check_Severe$Latitude.Degrees), breaks=break_pts_lat)

break_pts_lon<-seq(-180, 180, 10)
lon0<-hist(na.omit(Bleaching_Data_Reef_Check_None$Longitude.Degrees), breaks=break_pts_lon)
lon1<-hist(na.omit(Bleaching_Data_Reef_Check_Mild$Longitude.Degrees), breaks=break_pts_lon)
lon2<-hist(na.omit(Bleaching_Data_Reef_Check_Moderate$Longitude.Degrees), breaks=break_pts_lon)
lon3<-hist(na.omit(Bleaching_Data_Reef_Check_Severe$Longitude.Degrees), breaks=break_pts_lon)

#severity bleaching at Latitude
setwd(graphs_directory)
tiff(filename="Figure_S3_Final_color_blind.tif", res=400,width=2400,height=2400)
plot(lat0, xlab="Latitude", ylab="Number of surveys", main=NULL, col="#0D0887FF")
legend("topright", c("None", "Mild", "Moderate", "Severe"), fill=c("#0D0887FF", "#9C179EFF", "#ED7953FF", "#F0F921FF"), bty="n")
plot(lat1, col="#9C179EFF", add=T)
plot(lat2, col="#ED7953FF", add=T)
plot(lat3, col="#F0F921FF", add=T)
dev.off()

#Severity bleaching at Longitude
setwd(graphs_directory)
tiff(filename="Figure_S4_Final_color_blind.tif", res=400,width=3200,height=2400)
plot(lon0, xlab="Longitude", ylab="Number of surveys", main=NULL, col="#0D0887FF")
legend("topright", c("None", "Mild", "Moderate", "Severe"), fill=c("#0D0887FF", "#9C179EFF", "#ED7953FF", "#F0F921FF"), bty="n")
plot(lon1, col="#9C179EFF", add=T)
plot(lon2, col="#ED7953FF", add=T)
plot(lon3, col="#F0F921FF", add=T)
dev.off()

setwd(main_working_directory)
```

#Fig S3 with absolute value lat, for latitude
```{r plot severity bleaching per lat, lon, month, year}
#bleaching severity proportion per latitude, per longitude, per year, per month

break_pts_abs_lat<-seq(0, 40, 5)
lat_abs_0<-hist(na.omit(abs(Bleaching_Data_Reef_Check_None$Latitude.Degrees)), breaks=break_pts_abs_lat)
lat_abs_1<-hist(na.omit(abs(Bleaching_Data_Reef_Check_Mild$Latitude.Degrees)), breaks=break_pts_abs_lat)
lat_abs_2<-hist(na.omit(abs(Bleaching_Data_Reef_Check_Moderate$Latitude.Degrees)), breaks=break_pts_abs_lat)
lat_abs_3<-hist(na.omit(abs(Bleaching_Data_Reef_Check_Severe$Latitude.Degrees)), breaks=break_pts_abs_lat)


setwd(graphs_directory)
tiff(filename="bleaching_severity_frequency_histogram_abs_lat.tif", res=400,width=2400,height=2400)
plot(lat_abs_0, xlab="Latitude", ylab="Number of surveys", main=NULL, col="blue")
legend("topright", c("None", "Mild", "Moderate", "Severe"), fill=c("blue", "pink", "red", "black"), bty="n")
plot(lat_abs_1, col="pink", add=T)
plot(lat_abs_2, col="red", add=T)
plot(lat_abs_3, col="black", add=T)
dev.off()

plot(x=lat_abs_0$mids, y=lat_abs_0$counts, xlab="absolute value latitude", ylab="# surveys")
plot(x=lat_abs_0$mids, y=(lat_abs_0$counts-lat_abs_1$counts), xlab="absolute value latitude", ylab="# surveys with no bleaching")
plot(x=lat_abs_1$mids, y=lat_abs_1$counts-lat_abs_2$counts, xlab="absolute value latitude", ylab="# surveys with mild bleaching")
plot(x=lat_abs_2$mids, y=lat_abs_2$counts-lat_abs_3$counts, xlab="absolute value latitude", ylab="# surveys with moderate bleaching")
plot(x=lat_abs_3$mids, y=lat_abs_3$counts, xlab="absolute value latitude", ylab="# surveys with severe bleaching")
plot(x=lat_abs_3$mids, y=lat_abs_2$counts, xlab="absolute value latitude", ylab="# surveys with moderate or severe bleaching")

setwd(main_working_directory)
```

#figs S5-S15
```{r boxplots of variables at latitude}
plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$SSTA)
plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$SSTA_DHW)
plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$TSA)

plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$SSTA_Frequency)
plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$TSA_Frequency)
plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$Temperature_Kelvin)
plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$TSA_Maximum)
plot(Bleaching_Data_with_cortad_variables$Latitude.Degrees,Bleaching_Data_with_cortad_variables$Temperature_Kelvin_Standard_Deviation)

bin_breaks=c(-35, -30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35)

Bleaching_Data_with_cortad_variables$group=Bleaching_Data_with_cortad_variables$Latitude.Degrees
Bleaching_Data_with_cortad_variables$group=cut(Bleaching_Data_with_cortad_variables$group, breaks=bin_breaks, include.lowest=FALSE, right=TRUE, dig.lab=2, orderd_result=TRUE)

setwd(graphs_directory)
boxplot_function<-function(y, y_label){
  tiff(filename=paste(y_label, "_at_Latitude_boxplot.tif"), res=300,width=2000,height=2000)
  boxplot(y~group, data=Bleaching_Data_with_cortad_variables, las=2, xlab="", ylab=y_label)
  mtext(side=1, text="Latitude Range", line=4)
  dev.off()
}
boxplot_function(Bleaching_Data_with_cortad_variables$SSTA_Frequency, "SSTA_Frequency")
boxplot_function(Bleaching_Data_with_cortad_variables$SSTA_Frequency_Standard_Deviation, "SSTA_Frequency_stdev")
boxplot_function(Bleaching_Data_with_cortad_variables$TSA_Frequency, "TSA_Frequency")
boxplot_function(Bleaching_Data_with_cortad_variables$Temperature_Kelvin-273.15, "SST")
boxplot_function(Bleaching_Data_with_cortad_variables$TSA_Frequency, "TSA_Max")
boxplot_function(Bleaching_Data_with_cortad_variables$Temperature_Kelvin_Standard_Deviation, "SST_stdev")
boxplot_function(Bleaching_Data_with_cortad_variables$SSTA, "SSTA")
boxplot_function(Bleaching_Data_with_cortad_variables$SSTA_DHW, "SSTA_DHW")
boxplot_function(Bleaching_Data_with_cortad_variables$TSA, "TSA")
boxplot_function(Bleaching_Data_with_cortad_variables$TSA_Frequency_Standard_Deviation, "TSA_Freq_stdev")
boxplot_function(Bleaching_Data_with_cortad_variables$TSA_DHWMean, "TSA_DHW_Mean") #note* change ylim=c(0,4) to plot without the outliers
boxplot_function(Bleaching_Data_with_cortad_variables$TSA_DHW_Standard_Deviation, "TSA_DHW_stdev")
```

#fig S16
```{r covariate correlation plot}
MyVar <- c("Latitude.Degrees", "Year", "Depth", "Temperature_Kelvin", "Temperature_Kelvin_Standard_Deviation", "Temperature_Maximum", "SSTA", "SSTA_Maximum", "SSTA_Minimum", "SSTA_Frequency", "SSTA_Frequency_Standard_Deviation", "SSTA_DHW", "TSA_Frequency", "TSA_Frequency_Standard_Deviation", "TSA_DHW_Standard_Deviation", "ClimSST", "rate_of_SST_change")
#MyVar <- c("Latitude.Degrees.standardized", "Year.standardized", "Depth.standardized", "Temperature_Kelvin.standardized", "Temperature_Kelvin_Standard_Deviation.standardized", "Temperature_Maximum.standardized", "SSTA.standardized", "SSTA_Maximum.standardized", "SSTA_Minimum.standardized", "SSTA_Frequency.standardized", "SSTA_Frequency_Standard_Deviation.standardized", "SSTA_DHW.standardized", "TSA_Frequency.standardized", "TSA_Frequency_Standard_Deviation.standardized", "TSA_DHW_Standard_Deviation.standardized", "ClimSST.standardized", "rate_of_SST_change.standardized")

pairs(Bleaching_Data_with_cortad_variables[, MyVar], 
      lower.panel = panel.cor)

#variance-inflation factors; useful alongside pairplots and scatterplots
corvif(Bleaching_Data_with_cortad_variables[,MyVar])
M<-cor(Bleaching_Data_with_cortad_variables[,MyVar])
colnames(M) <- c("Latitude", "Year", "Depth", "SST", "SST_stdev", "SST_Max", "SSTA", "SSTA_Max", "SSTA_Min",  "SSTA_Freq", "SSTA_Freq_stdev", "SSTA_DHW", "TSA_Freq", "TSA_Freq_stdev", "TSA_DHW_stdev", "ClimSST", "Rate_of_SST_change")
rownames(M) <- c("Latitude", "Year", "Depth", "SST", "SST_stdev", "SST_Max", "SSTA", "SSTA_Max", "SSTA_Min",  "SSTA_Freq", "SSTA_Freq_stdev", "SSTA_DHW", "TSA_Freq", "TSA_Freq_stdev", "TSA_DHW_stdev", "ClimSST", "Rate_of_SST_change")
library(corrplot)
setwd(graphs_directory)
tiff(filename="Figure_S16_Final_color_blind.tif", res=300,width=2000,height=2000)
#corrplot(M, method="circle")
color_palette<-rev(plasma(100))
corrplot(M, method="circle", col=color_palette)
dev.off()
setwd(main_working_directory)
```

#figs S17, S19-S22
```{r plot variables in ecoregion}
for (i in 5:(dim(mean_variables_per_ecoregion)[2])){
setwd(graphs_directory)
mean_of_variable<-mean(mean_variables_per_ecoregion[,i])
sd_of_variable<-sd(mean_variables_per_ecoregion[,i])
standardized_value_of_variable<-(mean_variables_per_ecoregion[,i]-mean_of_variable)/sd_of_variable

D<-cbind(ERG=as.character(mean_variables_per_ecoregion$ERG), val=standardized_value_of_variable)
D<-as.data.frame(D)
D$val<-as.numeric(as.character(D$val))
nsteps <- 6;
step<-max(abs(max(na.omit(D$val))),abs(min(na.omit(D$val))))/nsteps
brks<-(-1)*max(abs(max(na.omit(D$val))),abs(min(na.omit(D$val))))+(0:(nsteps*2))*step
#cols <- c(rev(brewer.pal(9, "Blues")[4:9]), "gray75", (brewer.pal(9, "Reds")[4:9]))
cols <- c(plasma(100)[1], plasma(100)[5],plasma(100)[10],plasma(100)[15],plasma(100)[20],plasma(100)[35],plasma(100)[50],plasma(100)[65],plasma(100)[80],plasma(100)[85],plasma(100)[90],plasma(100)[95],plasma(100)[100])
grps<-(cut(D$val, brks, include.lowest=TRUE))
D$grp<-(cut(D$val, brks, include.lowest=TRUE))
ECO$grp<-NA
for(e in 1:length(ECO$grp)){
  if (ECO$ERG[e] %in%  D$ERG){
    ECO$grp[e]<-D$grp[(match(ECO$ERG, D$ERG))][e]
  }
}

tiff(paste("variable_mean_",names(mean_variables_per_ecoregion)[i],"_10_or_more_surveys_per_ecoregion_color_blind.tif", sep=""),res=300,width=2600,height=1000)
par(mai=c(.04,0.02,0.02,0.02))
library(maps)
plot.map("world", center=0 , ylim=c(-5500000,5500000), fill=T, col="grey90",border='grey70',mar=c(.1,.1,.1,.1)) #center is still 0
plot(ECO, ylim=c(-5500000,5500000), col=cols[unclass(ECO$grp)],add=T,lwd=.5)
axis(1,at=c(-10018754.17,3339584.724,16697920),lab=c('60°','180°','-60°'),las=1,cex.axis=.9,tcl=0.2,mgp=c(-1,-1.1,0))
axis(2, at=c(23*111319.4666666667,0,-23*111319.4666666667),labels=c('23°','0°','-23°'),las=3,cex.axis=.9,tcl=0.2,mgp=c(-1,-1.1,0),hadj=.4)
axis(3,at=c(-10018754.17,3339584.724,16697920),lab=c('','',''),las=1,tcl=0.2,mgp=c(-1,-1.3,0))
axis(4, at=c(23*111319.4666666667,0,-23*111319.4666666667),labels=c('','',''),las=2,tcl=0.2,mgp=c(-1,-0.6,0),hadj=0)
box()
windowsFonts(A=windowsFont('Arial Unicode MS'))
text(-7868896,-2922012,'Indian Ocean',cex=.9,family='A')
text(9438742,487176,'Pacific Ocean',cex=.9,family='A')
plotrix::color.legend(9684797.171+25e5,-28*111319.4666666667,15807371.62+25e5,-23.5*111319.4666666667,legend=c(round(min(brks), digits=1),round(max(brks), digits=1)),rect.col=cols,cex=.7)
north.arrow(x=(-16654136+111319.4*320), y=1615153*2, len=(111319.4*2), lab="N", cex=.8)
dev.off()
}
```

#Fig S18
```{r plot diversity}
setwd(graphs_directory)
D<-cbind(ERG=as.character(coral_diversity$Ecoregion), val=coral_diversity$SpeciesAccepted)
D<-as.data.frame(D)
D$val<-as.numeric(as.character(D$val))
nsteps <- 9;
step<-max(abs(max(na.omit(D$val))),abs(min(na.omit(D$val))))/nsteps
brks<-(0:(nsteps))*step
cols <- c(brewer.pal(9, "Reds"))
grps<-(cut(D$val, brks, include.lowest=TRUE))
D$grp<-(cut(D$val, brks, include.lowest=TRUE))
ECO$grp<-NA
for(e in 1:length(ECO$grp)){
  if (ECO$ERG[e] %in%  D$ERG){
    ECO$grp[e]<-D$grp[(match(ECO$ERG, D$ERG))][e]
  }
}

tiff(paste("Figure_S18_Final.tif", sep=""),res=300,width=2600,height=1000)
par(mai=c(.4,0.2,0.2,0.2))
library(maps)
plot.map("world", center=0 , ylim=c(-5500000,5500000), fill=T, col="grey90",border='grey70',mar=c(4,4,1,18))
plot(ECO, ylim=c(-5500000,5500000), col=cols[unclass(ECO$grp)],add=T,lwd=.5)
axis(1,at=c(-10018754.17,3339584.724,16697920),lab=c('60°','180°','-60°'),las=1,cex.axis=.7,tcl=0.2,mgp=c(-1,-1.1,0))
axis(2, at=c(23*111319.4666666667,0,-23*111319.4666666667),labels=c('23°','0°','-23°'),las=3,cex.axis=.7,tcl=0.2,mgp=c(-1,-1,0),hadj=.4)
axis(3,at=c(-10018754.17,3339584.724,16697920),lab=c('','',''),las=1,tcl=0.2,mgp=c(-1,-1.3,0))
axis(4, at=c(23*111319.4666666667,0,-23*111319.4666666667),labels=c('','',''),las=2,tcl=0.2,mgp=c(-1,-0.6,0),hadj=0)
box()
windowsFonts(A=windowsFont('Arial Unicode MS'))
text(-7868896,-2922012,'Indian Ocean',cex=.64,family='A')
text(9438742,487176,'Pacific Ocean',cex=.64,family='A')
plotrix::color.legend(9684797.171+25e5,-28*111319.4666666667,15807371.62+25e5,-23.5*111319.4666666667,legend=c(round(min(brks), digits=1),round(max(brks), digits=1)),rect.col=cols,cex=.7)
north.arrow(x=(-16654136+111319.4*320), y=1615153*2, len=(111319.4*2), lab="N", cex=.6)
dev.off()

```

#Fig S23
```{r bleaching presence and absence, early and late}
#2010 vs 2014 to 2017 event with 60 survey threshold
Bpresent <- subset(Bleaching_Data_with_cortad_variables, Average_bleaching >= 1, select=c(Temperature_Kelvin, Average_bleaching, Year)) 
Bpresent$Temp= Bpresent$Temperature_Kelvin - 273.15
Bpresentearly<-subset(Bpresent, Year==2010)  #
Bpresentlate<-subset(Bpresent, Year>=2014) #

BAbsent <-subset(Bleaching_Data_with_cortad_variables, Average_bleaching <1, select=c(Temperature_Kelvin, Average_bleaching, Year)) 
BAbsent$Temp= BAbsent$Temperature_Kelvin - 273.15
Babsentearly<-subset(BAbsent, Year==2010)
Babsentlate<-subset(BAbsent, Year>=2014)

presents<-hist(Bpresent$Temp, freq=T, breaks=c(seq(from=18, to=34, by=1)))
presentsearly<-hist(Bpresentearly$Temp, freq=T, breaks=c(seq(from=18, to=34, by=1)))
presentslate<-hist(Bpresentlate$Temp, freq=T, breaks=c(seq(from=18, to=34, by=1)))

absents<-hist(BAbsent$Temp, freq=T, breaks=c(seq(from=18, to=34, by=1)))
absentsearly<-hist(Babsentearly$Temp, freq=T, breaks=c(seq(from=18, to=34, by=1)))
absentslate<-hist(Babsentlate$Temp, freq=T, breaks=c(seq(from=18, to=34, by=1)))

xbins<-c(seq(from=18, to=34, by=1))
yvals<-presents$counts/(presents$counts+absents$counts)
yvalsearly<-presentsearly$counts/(presentsearly$counts+absentsearly$counts)
yvalslate<-presentslate$counts/(presentslate$counts+absentslate$counts)

xvals<-presents$mids

for (i in 1:length(yvals)){
  if(presentsearly$counts[i]+absentsearly$counts[i]<60){presentsearly$counts[i]<-0; absentsearly$counts[i]<-0}
  if(presentslate$counts[i]+absentslate$counts[i]<60){presentslate$counts[i]<-0; absentslate$counts[i]<-0}
}
yvalsearly<-presentsearly$counts/(presentsearly$counts+absentsearly$counts)
yvalslate<-presentslate$counts/(presentslate$counts+absentslate$counts)
setwd(graphs_directory)
tiff("Figure_S23_Final_color_blind.tif",res=300, width=1600, height=1600)
plot(x=xvals, y=yvalsearly, col="darkblue", xlab="°C", ylab="Proportion surveys with bleaching", ylim=c(0,.6))
legend("topleft", legend=c("2010 bleaching event", "2014-2017 bleaching event"), col=c("darkblue", "magenta"),pch=1, bty="o")
points(x=xvals, y=yvalslate, col="magenta")
dev.off()

```