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)

# Merge data

merged = merge(estacoes, rrs, by = 'BR_ID')

# change to numeric if necessary

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


# select data to plot

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))}

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


CHLA = merged %>% select('state', 'chla_ugL') %>% group_by(state) %>% summarise_each(fun = median.na)
ZSD = merged %>% select('state', 'secchi_m') %>% group_by(state) %>% summarise_each(fun = median.na)
ZSDzd = merged %>% select('state', 'secchi_m') %>% group_by(state) %>% summarise_each(fun = sd.na)

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

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


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


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

maximos = rrs_melt %>% group_by(region) %>% summarise(fun = max.na(max))


N$MAX = maximos$fun

rrs_melt$min[rrs_melt$min < 0] = 0

df = data.frame(region = N$state, 
                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()


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 = 4) + 
  
  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 = 8) +
  
  scale_y_continuous(name = expression(R[rs]~(sr^-1))) + 
  scale_x_continuous(name = "Wavelength (nm)", limits = c(400,900), 
                     breaks = seq(from = 400, to = 900, by = 100)) + 
  scale_fill_viridis_d(name = "region")+ 
  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_04.jpeg', 
       width = 30, height = 25, dpi = 200, units = 'in')