##############################################################################################################################
# Packages
##############################################################################################################################
library(data.table)

##############################################################################################################################
# 1. Data
##############################################################################################################################
# all points where CWCs occur have been sampled for environmental conditions in QGIS
# here, the values for environmental conditions at sites with CWCs are compared with the global mean

# set working directory to source file location

# Global Environmental Data
global_rawdata    <- read.csv("global_env.csv")
global            <- data.frame(global_rawdata, rep("global", 6183457))
colnames(global)  <- c("spec", colnames(global_rawdata[2:9]), "ID")


# Environmental Data at sites with Desmophyllum pertusum
despe_rawdata    <- read.csv("despe_env.csv")
despe              <- data.frame(despe_rawdata[, c(1, 5:12)], rep("despe", 8894))
colnames(despe)  <- c("spec", colnames(despe_rawdata[5:12]), "ID")

despe_depth1 <- data.frame(despe_rawdata$depth*-1, rep("despe", 8894))
colnames(despe_depth1) <- c("depth", "ID")


# Environmental Data at sites with Enallopsammia profunda
enapr_rawdata    <- read.csv("enapr_env.csv")
enapr            <- data.frame(enapr_rawdata[, c(1, 5:12)], rep("enapr", 505))
colnames(enapr)  <- c("spec", colnames(enapr_rawdata[5:12]), "ID")

enapr_depth1 <- data.frame(enapr_rawdata$depth*-1, rep("enapr", 505))
colnames(enapr_depth1) <- c("depth", "ID")


# Environmental Data at sites with Enallopsammia pusilla
enapu_rawdata    <- read.csv("enapu_env.csv")
enapu            <- data.frame(enapu_rawdata[, c(1, 5:12)], rep("enapu", 30))
colnames(enapu)  <- c("spec", colnames(enapu_rawdata[5:12]), "ID")

enapu_depth1 <- data.frame(enapu_rawdata$depth*-1, rep("enapu", 30))
colnames(enapu_depth1) <- c("depth", "ID")


# Environmental Data at sites with Enallopsammia rostrata
enaro_rawdata    <- read.csv("enaro_env.csv")
enaro            <- data.frame(enaro_rawdata[, c(1, 5:12)], rep("enaro", 1766))
colnames(enaro)  <- c("spec", colnames(enaro_rawdata[5:12]), "ID")

enaro_depth1 <- data.frame(enaro_rawdata$depth*-1, rep("enaro", 1766))
colnames(enaro_depth1) <- c("depth", "ID")


# Environmental Data at sites with Goniochorella dumosa
gondu_rawdata    <- read.csv("gondu_env.csv")
gondu            <- data.frame(gondu_rawdata[, c(1, 5:12)], rep("gondu", 718))
colnames(gondu)  <- c("spec", colnames(gondu_rawdata[5:12]), "ID")

gondu_depth1 <- data.frame(gondu_rawdata$depth*-1, rep("gondu", 718))
colnames(gondu_depth1) <- c("depth", "ID")

# Environmental Data at sites with Madrepora carolina
madca_rawdata    <- read.csv("madca_env.csv")
madca            <- data.frame(madca_rawdata[, c(1, 5:12)], rep("madca", 234))
colnames(madca)  <- c("spec", colnames(madca_rawdata[5:12]), "ID")

madca_depth1 <- data.frame(madca_rawdata$depth*-1, rep("madca", 234))
colnames(madca_depth1) <- c("depth", "ID")


# Environmental Data at sites with Madrepora oculata
madoc_rawdata    <- read.csv("madoc_env.csv")
madoc            <- data.frame(madoc_rawdata[, c(1, 5:12)], rep("madoc", 1524))
colnames(madoc)  <- c("spec", colnames(madoc_rawdata[5:12]), "ID")

madoc_depth1 <- data.frame(madoc_rawdata$depth*-1, rep("madoc", 1524))
colnames(madoc_depth1) <- c("depth", "ID")


# Environmental Data at sites with Solenosmilia vairiabilis
solva_rawdata    <- read.csv("solva_env.csv")
solva            <- data.frame(solva_rawdata[, c(1, 5:12)], rep("solva", 5574))
colnames(solva)  <- c("spec", colnames(solva_rawdata[5:12]), "ID")

solva_depth1 <- data.frame(solva_rawdata$depth*-1, rep("solva", 5574))
colnames(solva_depth1) <- c("depth", "ID")


##############################################################################################################################
# 2. Boxplots for Figure 3
##############################################################################################################################

cwc_spec       <- rbind(despe, enapr, enapu, enaro, gondu, madca, madoc, solva, global)
#cwc_spec_depth <- rbind(despe_depth1, enapr_depth1, enapu_depth1, enaro_depth1, gondu_depth1, 
#                        madca_depth1, madoc_depth1, solva_depth1, global_bath)

# Remove neg. values for current speed
cwc_spec_v <- cwc_spec[-which(cwc_spec$SAMPLE_v1<0),]
sum(cwc_spec_v$SAMPLE_v1<0)       # this is just to check if no more negative values

# Boxplots for Figure 3
plot.new()
par(mfrow = c(2,2), oma=c(0.5,0.45,0.1,0.1), mai = c(0.5,0.45,0.1,0.1), mgp=c(0,0.4,0), las=2)

boxplot(log(SAMPLE_PPmean1) ~ ID, data=cwc_spec, xlab="", ylab="", range=0)
title(ylab="log(PP_mean)", line=1.7, cex.lab=1.2, cex.lab=1)

boxplot(log(SAMPLE_PPrange1) ~ ID, data=cwc_spec, xlab="", ylab="", range=0)
title(ylab="log(PP_range)", line=1.7, cex.lab=1.2, cex.lab=1)

boxplot(log(SAMPLE_chl1) ~ ID, data=cwc_spec, xlab="", ylab="", range=0)    #this figure has not been used in the Artcile
title(ylab="log(chl)", line=1.7, cex.lab=1.2, cex.lab=1)

boxplot(log(SAMPLE_v1) ~ ID, data=cwc_spec_v, xlab="", ylab="", range=0)
title(ylab="log(v)", line=1.7, cex.lab=1.2, cex.lab=1)