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