require(dplyr)
require(ggplot2)
require(tidyr)
require(reshape2)
require(lemon)
require(plotly)
require(data.table)
require(geobr)
require(terra)
require(openxlsx)

estacoes = read.xlsx('Data/stations.xlsx', detectDates = T)

rrs = read.xlsx('Data/rrs.xlsx', detectDates = T)


rrs$Classe = fread('Outputs/SAM_OWT_rrs.csv')$Class

merged = merge(estacoes, rrs, by = 'BR_ID', no.dups = T)

for(i in 5:ncol(merged)) {
  
  merged[,i] = as.numeric(  merged[,i])
  
}

RRS = paste('Rrs_', 400:900, sep = '')

median.na = function(x) { return(median(x, na.rm = T))}
sd.na = function(x) { return(sd(x, na.rm = T))}
max.na = function(x) { return(max(x, na.rm = T))}
min.na = function(x) { return(min(x, na.rm = T))}

merged$estado = merged$Classe

merged = filter(merged, estado > 0)

res.mean = merged  %>% group_by(estado) %>% select(contains(RRS)) %>% summarise_each(fun = median.na)
res.sd =   merged  %>% group_by(estado) %>% select(contains(RRS)) %>% summarise_each(fun = sd.na)
max.sd =   merged  %>% group_by(estado) %>% select(contains(RRS)) %>% summarise_each(fun = max.na)
min.sd =   merged  %>% group_by(estado) %>% select(contains(RRS)) %>% summarise_each(fun = min.na)
N = merged %>% group_by(estado) %>% summarise(fun = n())

N$MAX = apply(max.sd[,-1], 1, max)
CHLA = merged %>% select('estado', 'chla_ugL') %>% group_by(estado) %>% summarise_each(fun = median.na)
ZSD = merged %>% select('estado', 'secchi_m') %>% group_by(estado) %>% summarise_each(fun = median.na)
ZSDzd = merged %>% select('estado', 'secchi_m') %>% group_by(estado) %>% summarise_each(fun = sd.na)

TSS = merged %>% select('estado', 'tss_mgL') %>% group_by(estado) %>% summarise_each(fun = median.na)
CDOM = merged %>% select('estado', 'acdom440') %>% group_by(estado) %>% summarise_each(fun = median.na)

rrs_melt.mean = res.mean %>% melt('estado')
rrs_melt.sd = res.sd %>% melt("estado")
rrs_melt.max = max.sd %>% melt("estado")
rrs_melt.min = min.sd %>% melt("estado")


rrs_melt = data.frame(media = rrs_melt.mean$value, 
                      region = rrs_melt.mean$estado,
                      variable = rrs_melt.mean$variable, 
                      sd = rrs_melt.sd$value,
                      max = rrs_melt.max$value,
                      min = rrs_melt.min$value)

df = data.frame(region = N$estado, 
                N = N$fun, 
                chla = round(CHLA$chla_ugL, 2),
                ZSD = round(ZSD$secchi_m, 2),
                TSS = round(TSS$tss_mgL, 2), 
                cdom = round(CDOM$acdom440, 2),
                ZSD.sd = round(ZSDzd$secchi_m, 2))

df$ZSD_plus_sd = paste(round(df$ZSD,2), '±', round(df$ZSD.sd,2))


rrs_melt$variable = gsub(rrs_melt$variable, pattern = "Rrs_", replacement = '') %>% as.numeric()

rrs_melt = filter(rrs_melt, region > 0)

rrs_melt$region = factor(x = rrs_melt$region, levels = c(1,2,3,4,5,6,7))

df = filter(df, region > 0)

df$region = factor(x = df$region, levels = c(1,2,3,4,5,6,7))



rrs_melt$region = gsub(rrs_melt$region, pattern = 1, replacement = 'OWT-01')
rrs_melt$region = gsub(rrs_melt$region, pattern = 2, replacement = 'OWT-02')
rrs_melt$region = gsub(rrs_melt$region, pattern = 3, replacement = 'OWT-03')
rrs_melt$region = gsub(rrs_melt$region, pattern = 4, replacement = 'OWT-04')
rrs_melt$region = gsub(rrs_melt$region, pattern = 5, replacement = 'OWT-05')
rrs_melt$region = gsub(rrs_melt$region, pattern = 6, replacement = 'OWT-06')
rrs_melt$region = gsub(rrs_melt$region, pattern = 7, replacement = 'OWT-07')

rrs_melt$region = factor(x = rrs_melt$region, levels = c('OWT-01',
                                                         'OWT-02',
                                                         'OWT-03',
                                                         'OWT-04',
                                                         'OWT-05',
                                                         'OWT-06',
                                                         'OWT-07'))

df$region = gsub(df$region, pattern = 1, replacement = 'OWT-01')
df$region = gsub(df$region, pattern = 2, replacement = 'OWT-02')
df$region = gsub(df$region, pattern = 3, replacement = 'OWT-03')
df$region = gsub(df$region, pattern = 4, replacement = 'OWT-04')
df$region = gsub(df$region, pattern = 5, replacement = 'OWT-05')
df$region = gsub(df$region, pattern = 6, replacement = 'OWT-06')
df$region = gsub(df$region, pattern = 7, replacement = 'OWT-07')

df$region = factor(x = df$region, levels = c('OWT-01',
                                             'OWT-02',
                                             'OWT-03',
                                             'OWT-04',
                                             'OWT-05',
                                             'OWT-06',
                                             'OWT-07'))


colorBlindBlack8  <- c("#0000FF", "#0072B2", "#56B4E9", "#009E73", 
                       "#F0E442", "#117733", "#D55E00")
size_axis = 20



plot_final = rrs_melt %>% ggplot(aes(x = variable)) + 
  geom_errorbar(aes(ymin = min, 
                    ymax = max, color = as.factor(region)), width=.2,
                position=position_dodge(0.05), alpha = 0.1) +
  
  geom_errorbar(aes(ymin = min, 
                    ymax = max+0.05, color = as.factor(region)), width=.2,
                position=position_dodge(0.05), alpha = 0) +
  
  geom_point(aes(y = media, color = as.factor(region)), size = 0.7) +
  
  
  facet_wrap(~(region), scales = 'free', ncol = 3) + 
  
  geom_text(x = 400, y = (N$MAX+0.05), hjust = 0, vjust = 1,
            aes(label = paste0("N = ", N,"\n", 
                               "Zsd = ", ZSD,"\n",
                               "Chl-a = ", chla,"\n",
                               "TSS = ", TSS,"\n",
                               "aCDOM(440) = ", cdom,"\n")), data = df, size = 7) +
  
  scale_y_continuous(name = expression(R[rs]~(sr^-1)), limits = c(NA,NA)) + 
  scale_x_continuous(name = "Wavelength (nm)", limits = c(400,900), 
                     breaks = seq(from = 400, to = 900, by = 100)) + 
  scale_fill_manual(values = colorBlindBlack8) +
  scale_color_manual(values = colorBlindBlack8) +
  theme_bw() + 
  theme(panel.grid.major = element_line(colour = "#d3d3d3"),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        text=element_text(family = "Tahoma"),
        axis.title = element_text( size = 35),
        axis.text.x = element_text(colour="black", size = size_axis),
        axis.text.y = element_text(colour="black", size = size_axis),
        axis.line = element_line(size=2, colour = "black"),
        strip.text = element_text(size = 40)) +
  #legend.title = element_text(size=20), #change legend title font size
  #legend.text = element_text(size=20),
  #legend.key.size = unit(2, 'cm')
  #change legend text font size) +
  theme(plot.margin = unit(c(4,4,4,4), "lines"))  + 
  guides(color = guide_legend(override.aes = list(size=15))) +
  theme(legend.position="none")


ggsave(plot = plot_final, filename = 'Outputs/Figures/Figure_05.jpeg', 
       width = 20, height = 17, dpi = 200, units = 'in')