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
estacoes$OWT_Pahlevan_class = rrs$Classe


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

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

  df = data.frame(owt = merged$Classe, state = merged$state) %>% 
              mutate(owt_estado = paste(owt, state)) %>% 
              group_by(owt_estado) %>% 
              summarise(fun = length(owt_estado)) %>% separate(owt_estado, into = c('OWT', "Estado"))

df = filter(df, as.numeric(OWT) > 0)

df$owt_estado = paste(df$Estado, df$OWT)
size_axis = 30

plot_final = df %>%
  ggplot(aes(y = fun, x = OWT, color = 'black', fill = (OWT))) + 
  geom_bar(stat = "identity", color = 'black', linewidth = 2) +#, position = "dodge") + # 'position = "dodge"' para barras separadas, se necessário
  facet_wrap(~Estado, scales = 'free_y', ncol = 3) + 
  labs(y = 'Frequency', x = 'OWT') +
  scale_fill_manual(values =  c("#0000FF", "#0072B2", "#56B4E9", "#009E73", 
                                                    "#F0E442", "#117733", "#D55E00")) +
  #scale_x_discrete(breaks = 1:7, labels = 1:7, expand = c(0, 0)) + # Expande o eixo X para garantir que apareçam todos os rótulos
  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)) + 
  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_06.jpeg', 
       width = 25, height = 30, dpi = 200, units = 'in')