Rivermouth Paper: Results display

1 Introduction and Methods

This interactive document summarizes the relationships between the size of a buffer radius around river mouths and metrics of coastal water clarity KdPAR and Kd490 (mean, median, variance, 50-10-90-95th percentiles, % of valid observations), based on the 2013-2022 period particularly looking at the effect of the river mouth classes. We suggest a new grouping of classes too.

We run Kruskall-Wallis on pair-wise comparisons (Wilcoxon rank sum tests with continuity correction) to look for significant differences in distribution of Kd490 and KdPAR values between classes and for each buffer sizes.

Data used:

• River names and river mouths location: River Environmental Classification (REC2) + Shane Orchard polish, January 2025 version.

• Coastal water clarity:

  • Downwelling light attenuation of PAR (400-700nm) coefficient (/m) (Gall et al. 2022), version 5, July 2002 to August 2023. https://gis.niwa.co.nz/portal/apps/experiencebuilder/template/?id=9794f29cd417493894df99d422c30ec2.

  • Downwelling light attenuation at 490nm coefficient (/m) OC-CCI product, 4km resolution, https://www.oceancolour.org/thredds/ncss/grid/CCI_ALL-v6.0-DAILY/dataset.html, using QAA to calculate the IOPs, then Lee (2005) to get Kd490 from IOPs.

Table 1: River mouth hydrosystem original classes (left) and suggested new grouping (right)

Classes
Waituna-type lagoon	Open Coast
Hāpua-type lagoon	Open Coast
Beach streams	Open Coast
Freshwater river mouth	Open Coast
Tidal river mouth	Open Coast
Tidal lagoon	Open Coast
Shallow drowned valley	Enclosed Coast
Deep drowned valley	Enclosed Coast
Fjord	Enclosed Coast

2 Results

library(tidyverse)
library(sf)
library(plotly)
library(viridis)
library(knitr)
library(broom)
library(DT)
library(ggpubr)
library(trend)

2.1 Figure 4 Kd490, KdPAR boxplot, all radii

kd_metrics <- read_csv('KdPAR_SCENZ_v05_Kd490_OCCCI_10yrMetrics_From_1D_AllBuffers_Rivermouths_REC2_class2.csv') %>% 
      select(-CLASS,Interface)


rivermouth <- read_csv('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Data/REC2rivermouths_Order5-8_Reg_BioR_EcoD_HSclass_named.csv') %>% 
  mutate(NAME_name2 = paste0(NAME,'_',name_2)) %>% 
  select(NAME_name2,CLASS,BioR__BioR,CATCHMT_AR,REGC2017_N,StreamOrde,RID) %>% 
  mutate(CLASS2 = case_when(
    CLASS == "2. Waituna-type lagoon" | CLASS == "3. Hāpua-type lagoon" | CLASS == "4. Beach stream" | CLASS == "5. Freshwater river mouth" | CLASS == "6. Tidal river mouth" | CLASS == "7. Tidal lagoon" ~ 'Open Coast',   
    CLASS == "8. Shallow drowned valley" | CLASS == "9. Deep drowned valley" | CLASS == "10. Fjord" ~ 'Enclosed Coast',
    TRUE ~ 'NA')) %>% 
    mutate(CLASS = case_when(
    CLASS == "2. Waituna-type lagoon" ~ "Waituna-type lagoon",
    CLASS == "3. Hāpua-type lagoon" ~ "Hāpua-type lagoon",
    CLASS == "4. Beach stream" ~ "Beach stream",
    CLASS == "5. Freshwater river mouth" ~ "Freshwater river mouth",
    CLASS == "6. Tidal river mouth" ~ "Tidal river mouth",
    CLASS == "7. Tidal lagoon" ~ 'Tidal lagoon',   
    CLASS == "8. Shallow drowned valley" ~ "Shallow drowned valley",
    CLASS == "9. Deep drowned valley" ~ "Deep drowned valley",
    CLASS == "10. Fjord" ~ 'Fjord',
    TRUE ~ 'NA')) %>% 
    mutate(CLASS = factor(CLASS, levels = c("Waituna-type lagoon", 
                                          "Hāpua-type lagoon", 
                                          "Beach stream",
                                          "Freshwater river mouth", 
                                          "Tidal river mouth", 
                                          "Tidal lagoon",
                                          "Shallow drowned valley",
                                          "Deep drowned valley",
                                          "Fjord")))

kd_metrics <- inner_join(kd_metrics,rivermouth)

kd_variance <- kd_metrics %>% group_by(Buffer_Size) %>% 
  filter(Metrics == 'mean') %>% 
  mutate(variance = var(KdPAR,na.rm=T),
         Buffer_Size = Buffer_Size/1000) %>% 
  filter(!is.na(KdPAR))

ratio <- 15

p <- ggplot(data=kd_variance) +
  geom_boxplot(aes(x=factor(Buffer_Size),y=KdPAR)) +
  geom_path(aes(x= factor(Buffer_Size), y=variance*ratio), colour='red',group=1) +
  geom_point(aes(x= factor(Buffer_Size), y=variance*ratio), colour='red',group=1) +
  scale_y_continuous(
    name = "KdPAR (/m)",
    # Add a second axis and specify its features
    sec.axis = sec_axis(~ ./ratio, name="Variance (/m)^2")) +
  scale_fill_viridis_d() +
  xlab('Radius (km)') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.y.right = element_text(color = "red"),
        axis.line.y.right = element_line(color = "red"),
        axis.title.y.right = element_text(color = "red"))

kd_variance2 <- kd_metrics %>% group_by(Buffer_Size) %>% 
  filter(Metrics == 'mean') %>% 
  mutate(variance = var(Kd490,na.rm=T),
         Buffer_Size = Buffer_Size/1000) %>% 
  filter(!is.na(Kd490))

ratio <- 60

p2 <- ggplot(data=kd_variance2) +
  geom_boxplot(aes(x=factor(Buffer_Size),y=Kd490)) +
  geom_path(aes(x= factor(Buffer_Size), y=variance*ratio), colour='red',group=1) +
  geom_point(aes(x= factor(Buffer_Size), y=variance*ratio), colour='red',group=1) +
  scale_y_continuous(
    name = "Kd490 (/m)",
    # Add a second axis and specify its features
    sec.axis = sec_axis(~ ./ratio, name="Variance (/m)^2")) +
  scale_fill_viridis_d() +
  xlab('Radius (km)') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.y.right = element_text(color = "red"),
        axis.line.y.right = element_line(color = "red"),
        axis.title.y.right = element_text(color = "red"))

p_obs <- kd_metrics %>% 
  filter(Metrics == 'valid_days_pc') %>% 
  filter(!is.na(KdPAR)) %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  ggplot(data=.,aes(x=factor(Buffer_Size),y=KdPAR)) +
  geom_boxplot() +
  scale_fill_viridis_d() +
  #ggtitle('Percentage of valid observations of KdPAR (2013-2022)') +
  ylab('Percent of valid observations (%)') + xlab('Radius (km)') +
  theme_light() +
  theme(legend.position = 'bottom')

p2_obs <- kd_metrics %>% 
  filter(Metrics == 'valid_days_pc') %>% 
  filter(!is.na(Kd490)) %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  ggplot(data=.,aes(x=factor(Buffer_Size),y=Kd490)) +
  geom_boxplot() +
  scale_fill_viridis_d() +
  #ggtitle('Percentage of valid observations of Kd490 (2013-2022)') +
  ylab('Percent of valid observations (%)') + xlab('Radius (km)') +
  theme_light() +
  theme(legend.position = 'bottom')

figure1 <- ggarrange(p2,
                     p,
                     p2_obs,
                     p_obs,
                     labels = c("A", "B","C","D"),
                     ncol = 2, nrow = 2, align = 'v')
figure1

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure4_Kd_vs_Buffer_ValidObs.png',width = 8, height = 10, dpi=600, figure1)

2.2 Table 1 KdPAR and Kd490 summary per radius

table1_kdPAR <- kd_metrics %>% 
  filter(!is.na(KdPAR)) %>% 
  filter(Metrics == 'mean') %>% 
  group_by(Buffer_Size) %>% 
  summarise(mean = mean(KdPAR,na.rm=T),
         sd = sd(KdPAR,na.rm=T),
         se = sd/sqrt(length(KdPAR))) %>% 
    mutate(Buffer_Size = Buffer_Size/1000)

DT::datatable(table1_kdPAR) %>% 
  formatRound(columns = c('mean','sd','se'), digits=3)

2.2.1 KdPAR Pairwise comparison

aa <-  kd_metrics %>% 
  filter(!is.na(KdPAR)) %>% 
  filter(Metrics == 'mean') %>% 
  select(KdPAR,Buffer_Size) %>% 
  mutate(Buffer_Size = Buffer_Size/1000)


pw <- pairwise.wilcox.test(aa$KdPAR,aa$Buffer_Size)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('1',
                    '2',
                    '3',
                    '4',
                    '5',
                    '10',
                    '15'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatRound(c(2:8), 5) %>% 
  formatStyle(columns = c(1:8), 'text-align' = 'center')

table1_kd490 <- kd_metrics %>% 
  filter(!is.na(KdPAR)) %>% 
  filter(Metrics == 'mean') %>% 
  group_by(Buffer_Size) %>% 
  summarise(mean = mean(Kd490,na.rm=T),
         sd = sd(Kd490,na.rm=T),
         se = sd/sqrt(length(KdPAR))) %>% 
    mutate(Buffer_Size = Buffer_Size/1000) %>% 
  na.omit()

DT::datatable(table1_kd490) %>% 
  formatRound(columns = c('mean','sd','se'), digits=3)

2.2.2 Kd490 Pairwise comparison

aa2 <-  kd_metrics %>% 
  filter(!is.na(Kd490)) %>% 
  filter(Metrics == 'mean') %>% 
  select(Kd490,Buffer_Size) %>% 
  mutate(Buffer_Size = Buffer_Size/1000)


pw2 <- pairwise.wilcox.test(aa2$Kd490,aa2$Buffer_Size)

pw2 %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5',
                    '10',
                    '15',
                    '20'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatRound(c(2:5), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.3 Figure 5 Kd490 and KdPAR, boxplot, 10km buffer, coastal hydrosystem types (10)

p <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  ggplot(data=.) +
  geom_boxplot(aes(x=factor(CLASS),y=KdPAR)) +
  scale_fill_viridis_d() +
  xlab('') +
  ylab('KdPAR (/m)') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))


p2 <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  ggplot(data=.) +
  geom_boxplot(aes(x=factor(CLASS),y=Kd490)) +
  scale_fill_viridis_d() +
  xlab('') +
  ylab('Kd490 (/m)') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

figure2 <- ggarrange(p2,
                     p,
                     labels = c("A", "B"),
                     ncol = 2, nrow = 1)
figure2

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure5_Kd_vs_Hydrosystems.png',width = 6, height = 4, dpi=600, figure2)#8,6

2.4 Table 2 KdPAR and Kd490 summary per class, 10km radius

table2_kdPAR <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  filter(!is.na(KdPAR)) %>%
  group_by(CLASS) %>% 
  summarise(mean = mean(KdPAR,na.rm=T),
         sd = sd(KdPAR,na.rm=T),
         se = sd/sqrt(length(KdPAR))) 

DT::datatable(table2_kdPAR) %>% 
  formatRound(columns = c('mean','sd','se'), digits=3)

2.4.1 KdPAR Pairwise comparison

aa <- kd_metrics %>%   
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  select(KdPAR,CLASS) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$CLASS)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('Waituna-type lagoon',
                    'Hāpua-type lagoon',
                    'Beach stream',
                    'Freshwater river mouth',
                    'Tidal river mouth',
                    'Tidal lagoon',
                    'Shallow drowned valley',
                    'Deep drowned valley'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatRound(c(2:9), 5) %>% 
  formatStyle(columns = c(1:9), 'text-align' = 'center')

table2_kd490 <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  filter(!is.na(Kd490)) %>%
  group_by(CLASS) %>% 
  summarise(mean = mean(Kd490,na.rm=T),
         sd = sd(Kd490,na.rm=T),
         se = sd/sqrt(length(Kd490))) 

DT::datatable(table2_kd490) %>% 
  formatRound(columns = c('mean','sd','se'), digits=3)

2.4.2 Kd490 Pairwise comparison

aa2 <- kd_metrics %>%   
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  select(Kd490,CLASS) %>% na.omit()

pw2 <- pairwise.wilcox.test(aa2$Kd490,aa2$CLASS)

pw2 %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('Waituna-type lagoon',
                    'Hāpua-type lagoon',
                    'Beach stream',
                    'Freshwater river mouth',
                    'Tidal river mouth',
                    'Tidal lagoon',
                    'Shallow drowned valley'),
                  target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatRound(c(2:8), 5) %>% 
  formatStyle(columns = c(1:8), 'text-align' = 'center')

2.5 Figure 6 Kd490 and KdPAR Boxplot, 10yr mean, 10km buffer, NEW grouping

p <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  ggplot(data=.) +
  geom_boxplot(aes(x=factor(CLASS2),y=KdPAR)) +
  scale_fill_viridis_d() +
  xlab('') +
  ylab('KdPAR (/m)') +
  theme_light() +
  theme(legend.position = 'bottom')

p2 <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  ggplot(data=.) +
  geom_boxplot(aes(x=factor(CLASS2),y=Kd490)) +
  scale_fill_viridis_d() +
  xlab('') +
  ylab('Kd490 (/m)') +
  theme_light() +
  theme(legend.position = 'bottom')#axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)

figure3 <- ggarrange(p2,
                     p,
                     labels = c("A", "B"),
                     ncol = 2, nrow = 1)
figure3

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure6_Kd_vs_Geometry.png',width = 6, height = 4, dpi=600, figure3)

2.5.1 KdPAR Pairwise comparison

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  select(KdPAR,CLASS2) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$CLASS2)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('Enclosed Coast'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:2), 'text-align' = 'center')

2.5.2 Kd490 Pairwise comparison

aa2 <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'mean') %>% 
  select(Kd490,CLASS2) 

pw2 <- pairwise.wilcox.test(aa2$Kd490,aa2$CLASS2)

pw2 %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('Enclosed Coast'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatRound(c(2), 3) %>% 
  formatStyle(columns = c(1:2), 'text-align' = 'center')

2.6 Figure 7 Does stream order influence KdPAR? Analysis of 6 metrics

2.6.0.1 Main figure of case study

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & CLASS2 == "Open Coast" & Metrics %in% c('percentile10','median','percentile90')) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
    mutate(Metrics = factor(Metrics,levels = c('percentile5','10th percentile','Median','90th percentile','percentile95'))) %>% 
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~Metrics,ncol=3, scales = 'free') +
  scale_fill_viridis_d() +
  ylab('KdPAR (/m)') + xlab('Stream Order') +
  ggtitle("Open Coast") +
  theme_light() +
  theme(legend.position = 'bottom')
p

p2 <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & CLASS2 == "Enclosed Coast" & Metrics %in% c('percentile10','median','percentile90')) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
    mutate(Metrics = factor(Metrics,levels = c('percentile5','10th percentile','Median','90th percentile','percentile95'))) %>% 
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~Metrics,ncol=3, scales = 'free') +
  scale_fill_viridis_d() +
  ylab('KdPAR (/m)') + xlab('Stream Order') +
  ggtitle("Enclosed Coast") +
  theme_light() +
  theme(legend.position = 'bottom')
p2

figure7 <- ggarrange(p,
                     p2,
                     labels = c("A", "B"),
                     ncol = 1, nrow = 2)
figure7

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure7_Kd_vs_StreamOrder.png',width = 6, height = 4, dpi=600, figure7)

2.6.0.2 Percentile 5

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile5') %>%
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 5th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

2.6.0.2.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile5' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.2.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile5' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.3 Percentile 10

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile10') %>% 
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 10th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

2.6.0.3.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile10' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.3.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile10' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.4 Mean

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'mean') %>% 
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) Mean (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

2.6.0.4.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'mean' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.4.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'mean' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.5 Median

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'median') %>% 
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) Median (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

2.6.0.5.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'median' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.5.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'median' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.6 Percentile 90

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile90') %>% 
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 90th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

2.6.0.6.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile90' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.6.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile90' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.7 Percentile 95

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile95') %>% 
  ggplot(data=.,aes(x=factor(StreamOrde),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 95th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

2.6.0.7.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile95' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.6.0.7.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
           Metrics == 'percentile95' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,StreamOrde) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$StreamOrde)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(c('5', '6','7'),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:5), 'text-align' = 'center')

2.7 Figure 8 Does Bioregions influence KdPAR? Analysis of 6 metrics

2.7.0.1 Main figure of case study

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & CLASS2 == "Open Coast" & Metrics %in% c('percentile10','median','percentile90')) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
    mutate(Metrics = factor(Metrics,levels = c('percentile5','10th percentile','Median','90th percentile','percentile95'))) %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~Metrics,ncol=3,scales = 'free') +
  scale_fill_viridis_d() +
  ylab('KdPAR (/m)') + xlab('Bioregions') +
  ggtitle("Open Coast") +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

p2 <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & CLASS2 == "Enclosed Coast" & Metrics %in% c('percentile10','median','percentile90')) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
    mutate(Metrics = factor(Metrics,levels = c('percentile5','10th percentile','Median','90th percentile','percentile95'))) %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~Metrics,ncol=3,scales = 'free') +
  scale_fill_viridis_d() +
  ylab('KdPAR (/m)') + xlab('Bioregions') +
  ggtitle("Enclosed Coast") +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p2

figure8 <- ggarrange(p,
                     p2,
                     labels = c("A", "B"),
                     ncol = 1, nrow = 2)
figure8

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure8_Kd_vs_Bioregions.png',width = 12, height = 8, dpi=600, figure8)

2.7.0.2 Percentile 5

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile5') %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 5th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.7.0.2.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile5' &
         CLASS2 == 'Open Coast') %>% 
  #mutate(CLASS2_BioR = paste0(CLASS2,"_",BioR__BioR)) %>% 
  select(KdPAR,BioR__BioR) 

#pw <- pairwise.wilcox.test(aa$KdPAR,aa$CLASS2_BioR)
pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.2.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile5' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  #formatRound(c(2), 5) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.3 Percentile 10

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile10') %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 10th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.7.0.3.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile10' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.3.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile10' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.4 Mean

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'mean') %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) Mean (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.7.0.4.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'mean' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.4.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'mean' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.5 Median

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'median') %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) Median (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.7.0.5.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'median' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.5.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'median' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.6 Percentile 90

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile90') %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 90th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.7.0.6.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile90' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.6.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile90' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.7 Percentile 95

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile95') %>% 
  ggplot(data=.,aes(x=factor(BioR__BioR),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 95th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.7.0.7.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile95' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.7.0.7.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile95' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,BioR__BioR) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$BioR__BioR)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8 Figure XX Does Council regions influence KdPAR? Analysis of 6 metrics

2.8.0.1 Percentile 5

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile5') %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 5th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.8.0.1.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile5' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.1.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile5' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.2 Percentile 10

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile10') %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 10th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.8.0.2.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile10' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.2.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile10' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.3 Mean

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'mean') %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) mean (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.8.0.3.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'mean' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.3.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'mean' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.4 Median

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'median') %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) median (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.8.0.4.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'median' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.4.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'median' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.5 Percentile 90

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile90') %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 90th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.8.0.5.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile90' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.5.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile90' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.6 Percentile 95

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile95') %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  scale_fill_viridis_d() +
  ylab('KdPAR 10yr (2013-2022) 95th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

2.8.0.6.1 Open Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile95' &
         CLASS2 == 'Open Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.8.0.6.2 Enclosed Coast

aa <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & 
         Metrics == 'percentile95' &
         CLASS2 == 'Enclosed Coast') %>% 
  select(KdPAR,REGC2017_N) 

pw <- pairwise.wilcox.test(aa$KdPAR,aa$REGC2017_N)

pw %>% 
  tidy() %>% 
  pivot_wider(names_from = c(group2),values_from = p.value) %>% 
  DT::datatable() %>% 
  DT::formatStyle(columns = c(1:16),target = 'cell',fontWeight = styleInterval(.05, c('bold', 'normal'))) %>% 
  formatStyle(columns = c(1:16), 'text-align' = 'center')

2.9 Figure 9 Time series of National and river with highest KdPAR

2.9.1 Time series national scale

kd_ts_10km <- read_csv('KdPAR_SCENZ_1D_v05_Buffer_10000m_RiverMouths_REC2.csv')

kd_ts_10km <- inner_join(kd_ts_10km,rivermouth %>% select(RID,NAME_name2,CLASS2))

kd_ts_10km_yearly_sum <- kd_ts_10km %>% 
  mutate(year = year(Date)) %>% 
  group_by(year,CLASS2) %>% 
  summarise(KdPAR_median = median(KdPAR,na.rm=T),
            KdPAR_10p = quantile(KdPAR,na.rm=T,0.1),
            KdPAR_90p = quantile(KdPAR,na.rm=T,0.9))
  
p <- kd_ts_10km_yearly_sum %>% 
  filter(CLASS2 == 'Open Coast') %>% 
  ggplot(data=.,aes(x=year,y=KdPAR_median)) +
  geom_point() +
  geom_line() + 
  ylab('KdPAR (/m)') + xlab('Year') +
  geom_smooth(method = 'lm') +
  ggtitle("Open Coast") +
  #facet_wrap(~CLASS2,nrow=2,scales = 'free') +
  theme_light()
p

p2 <- kd_ts_10km_yearly_sum %>% 
  filter(CLASS2 == 'Enclosed Coast') %>% 
  ggplot(data=.,aes(x=year,y=KdPAR_median)) +
  geom_point() +
  geom_line() + 
  ylab('KdPAR (/m)') + xlab('Year') +
  geom_smooth(method = 'lm') +
  ggtitle("Enclosed Coast") +
  #facet_wrap(~CLASS2,nrow=2,scales = 'free') +
  theme_light()
p2

figure9 <- ggarrange(p,
                     p2,
                     labels = c("A", "B"),
                     ncol = 1, nrow = 2)
#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure9_ts_national.png',width = 6, height = 4, dpi=600, figure9)

p <- kd_ts_10km_yearly_sum %>% 
  filter(CLASS2 == 'Open Coast') %>% 
  pivot_longer(-c(year,CLASS2)) %>% 
  mutate(name = case_when(name == 'KdPAR_10p' ~  '10th percentile',
                             name == 'KdPAR_median' ~  'Median',
                             name == 'KdPAR_90p' ~  '90th percentile')) %>%
  mutate(name = factor(name,levels = c('10th percentile','Median','90th percentile'))) %>% 
  ggplot(data=.,aes(x=year,y=value)) +
  geom_point() +
  geom_line() + 
  ylab('KdPAR (/m)') + xlab('Year') +
  geom_smooth(method = 'lm') +
  ggtitle("Open Coast") +
  facet_wrap(~name,nrow=1,scales = 'free') +
  theme_light()
p

`geom_smooth()` using formula = 'y ~ x'

p2 <- kd_ts_10km_yearly_sum %>% 
  filter(CLASS2 == 'Enclosed Coast') %>% 
  pivot_longer(-c(year,CLASS2)) %>% 
  mutate(name = case_when(name == 'KdPAR_10p' ~  '10th percentile',
                             name == 'KdPAR_median' ~  'Median',
                             name == 'KdPAR_90p' ~  '90th percentile')) %>%
  mutate(name = factor(name,levels = c('10th percentile','Median','90th percentile'))) %>% 
  ggplot(data=.,aes(x=year,y=value)) +
  geom_point() +
  geom_line() + 
  ylab('KdPAR (/m)') + xlab('Year') +
  geom_smooth(method = 'lm') +
  ggtitle("Enclosed Coast") +
  facet_wrap(~name,nrow=1,scales = 'free') +
  theme_light()
p2

`geom_smooth()` using formula = 'y ~ x'

figure9b <- ggarrange(p,
                     p2,
                     labels = c("A", "B"),
                     ncol = 1, nrow = 2)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure9_ts_national_quantile.png',width = 6, height = 4, dpi=600, figure9b)

2.9.1.1 Trends

trends <- kd_ts_10km_yearly_sum %>% 
  pivot_longer(-c(year,CLASS2)) %>% 
  group_by(CLASS2,name) %>% 
  nest() %>% 
    mutate(ts_out = purrr::map(data, ~ts(.x$value,start=2002,end=2023,frequency = 1))) %>% 
  mutate(sens = purrr::map(ts_out, ~sens.slope(.x, conf.level = 0.95))) %>% 
  mutate(Sens_Slope = as.numeric(unlist(sens)[1]),P_Value =as.numeric(unlist(sens)[3])) %>% 
  select(CLASS2,name, Sens_Slope, P_Value)
  
write_csv(trends,'Trends_Figure9_ts_national_quantile.csv')

trends %>% 
  DT::datatable() 

2.9.2 Rivermouths ranked by highest KdPAR values (mean 2002-2023 period)

## Rankings of most dirty rivermouth
kd_ts_10km_summary <- kd_ts_10km %>% 
  group_by(RID,NAME,name_2,CLASS2) %>% 
  summarise(KdPAR_mean = round(mean(KdPAR,na.rm=T),3),
            KdPAR_sd = round(sd(KdPAR,na.rm=T),3)) %>% 
  arrange(desc(KdPAR_mean))

`summarise()` has grouped output by 'RID', 'NAME', 'name_2'. You can override
using the `.groups` argument.

#write_csv(kd_ts_10km_summary,'KdPAR_2002_2023_mean_sd_Rivermouth_ranked.csv')

kd_ts_10km_summary %>% 
  DT::datatable() 

2.9.3 Timeseries of rivermouth with highest KdPAR values (Kaipara Harbour System)

kd_ts_10km_yearly_sum_kai <- kd_ts_10km %>% 
  filter(RID == 12) %>% 
  mutate(year = year(Date)) %>% 
  group_by(year,CLASS2) %>% 
  summarise(KdPAR_median = median(KdPAR,na.rm=T),
            KdPAR_10p = quantile(KdPAR,na.rm=T,0.1),
            KdPAR_90p = quantile(KdPAR,na.rm=T,0.9))

`summarise()` has grouped output by 'year'. You can override using the
`.groups` argument.

p <- ggplot(data=kd_ts_10km_yearly_sum_kai,aes(x=year,y=KdPAR_median)) +
  geom_point() +
  geom_line() + 
  ylab('KdPAR (/m)') + xlab('Year') +
  geom_smooth(method = 'lm') +
  facet_wrap(~CLASS2,nrow=2,scales = 'free') +
  theme_light()
p

`geom_smooth()` using formula = 'y ~ x'

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure9_ts_kaipara.png',width = 5, height = 4, dpi=600, p)

2.10 Figure 10 KdPAR mean for all rivers

kd_mean_10km <- kd_ts_10km %>% group_by(CLASS2,RID,NAME,name_2,NAME_name2) %>% 
  summarise(KdPAR_mean = mean(KdPAR,na.rm=T)) %>% 
  arrange(desc(KdPAR_mean))

`summarise()` has grouped output by 'CLASS2', 'RID', 'NAME', 'name_2'. You can
override using the `.groups` argument.

p <- kd_mean_10km %>%
  arrange(desc(KdPAR_mean)) %>% 
  ggplot(data=.,aes(x=fct_inorder(NAME_name2) ,y=KdPAR_mean,colour=CLASS2)) +
  geom_segment(aes(xend=NAME_name2 ,y = 0, yend=KdPAR_mean,colour=CLASS2),size=.5) +
  scale_colour_manual(values=c('blue','red')) +
  geom_point(size=2) +
  xlab('Rivers (n = 226)') + ylab('KdPAR mean (/m)') +
  theme_classic() +
  theme(legend.position = 'bottom',
        legend.title = element_blank(),
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank())

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

p

p2 <- kd_mean_10km[1:20,] %>%
  ggplot(data=.,aes(x=fct_inorder(name_2) ,y=KdPAR_mean,colour=CLASS2)) +
  geom_segment(aes(xend= ,y = 0, yend=KdPAR_mean,colour=CLASS2)) +
  scale_colour_manual(values=c('blue','red')) +
  geom_point(size=2) +
  xlab('Rivers (top 20)') + ylab('KdPAR mean (/m)') +
  theme_light() +
  theme(legend.position = 'none',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p2

figure10 <- ggarrange(p,
                     p2,
                     widths = c(3, 2),
                     labels = c("A", "B"),
                     ncol = 2, nrow = 1)
figure10

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/Figure10_KdPAR_Mean_rivers_v2.png',width = 10, height = 6, dpi=600, figure10)

2.11 Figure XX Does Catchment area influence KdPAR? Analysis of 6 metrics

2.11.0.1 Percentile 5

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile5' & CATCHMT_AR>0) %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  ylab('KdPAR 10yr (2013-2022) 5th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

aa <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile5' & CATCHMT_AR>0) %>% 
  select(KdPAR, CATCHMT_AR, CLASS2) %>% 
  group_by(CLASS2) %>% 
  nest() %>% 
  mutate(model = map(data, ~broom::glance(lm(KdPAR~CATCHMT_AR, data=.x)))) %>% 
  select(-data) %>% 
  unnest(model)  

aa %>% DT::datatable()

2.11.0.2 Percentile 10

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile10' & CATCHMT_AR>0) %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  ylab('KdPAR 10yr (2013-2022) 10th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

aa <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile10' & CATCHMT_AR>0) %>% 
  select(KdPAR, CATCHMT_AR, CLASS2) %>% 
  group_by(CLASS2) %>% 
  nest() %>% 
  mutate(model = map(data, ~broom::glance(lm(KdPAR~CATCHMT_AR, data=.x)))) %>% 
  select(-data) %>% 
  unnest(model)  

aa %>% DT::datatable()

2.11.0.3 Mean

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'mean' & CATCHMT_AR>0) %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  ylab('KdPAR 10yr (2013-2022) Mean (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

aa <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'mean' & CATCHMT_AR>0) %>% 
  select(KdPAR, CATCHMT_AR, CLASS2) %>% 
  group_by(CLASS2) %>% 
  nest() %>% 
  mutate(model = map(data, ~broom::glance(lm(KdPAR~CATCHMT_AR, data=.x)))) %>% 
  select(-data) %>% 
  unnest(model)  

aa %>% DT::datatable()

2.11.0.4 Median

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'median' & CATCHMT_AR>0) %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  ylab('KdPAR 10yr (2013-2022) Median (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

aa <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'median' & CATCHMT_AR>0) %>% 
  select(KdPAR, CATCHMT_AR, CLASS2) %>% 
  group_by(CLASS2) %>% 
  nest() %>% 
  mutate(model = map(data, ~broom::glance(lm(KdPAR~CATCHMT_AR, data=.x)))) %>% 
  select(-data) %>% 
  unnest(model)  

aa %>% DT::datatable()

2.11.0.5 Percentile 90

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile90' & CATCHMT_AR>0) %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  ylab('KdPAR 10yr (2013-2022) 90th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

aa <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile90' & CATCHMT_AR>0) %>% 
  select(KdPAR, CATCHMT_AR, CLASS2) %>% 
  group_by(CLASS2) %>% 
  nest() %>% 
  mutate(model = map(data, ~broom::glance(lm(KdPAR~CATCHMT_AR, data=.x)))) %>% 
  select(-data) %>% 
  unnest(model)  

aa %>% DT::datatable()

2.11.0.6 Percentile 95

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile95' & CATCHMT_AR>0) %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  facet_wrap(~CLASS2,scales='free',ncol=2) +
  ylab('KdPAR 10yr (2013-2022) 95th percentile (/m)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom')
p

aa <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile95' & CATCHMT_AR>0) %>% 
  select(KdPAR, CATCHMT_AR, CLASS2) %>% 
  group_by(CLASS2) %>% 
  nest() %>% 
  mutate(model = map(data, ~broom::glance(lm(KdPAR~CATCHMT_AR, data=.x)))) %>% 
  select(-data) %>% 
  unnest(model)  

aa %>% DT::datatable()

3 Supplementary Information

3.1 Fig S1

figS1 <- kd_metrics %>% 
  filter(Metrics %in% c('median','percentile10','percentile90')) %>% 
  filter(!is.na(KdPAR)) %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
  mutate(Metrics = factor(Metrics,levels = c('10th percentile','Median','90th percentile'))) %>% 
  ggplot(data=.,aes(x=factor(Buffer_Size),y=KdPAR,fill=Metrics)) +
  geom_boxplot() +
  scale_fill_viridis_d() +
  facet_wrap(~CLASS) +
  #ggtitle('Percentage of valid observations of KdPAR (2013-2022)') +
  ylab('KdPAR (/m)') + xlab('Radius (km)') +
  theme_light() +
  theme(legend.position = 'bottom')
figS1

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/FigureS1_KdPAR_boxplot_CLASS_vs_Radii.png',width = 8, height = 6, dpi=600, figS1)

3.2 Fig S2

figS2 <- kd_metrics %>% 
  filter(Metrics %in% c('median','percentile10','percentile90')) %>% 
  filter(!is.na(KdPAR)) %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
  mutate(Metrics = factor(Metrics,levels = c('10th percentile','Median','90th percentile'))) %>% 
  ggplot(data=.,aes(x=factor(Buffer_Size),y=KdPAR,fill=Metrics)) +
  geom_boxplot() +
  scale_fill_viridis_d() +
  facet_wrap(~CLASS2) +
  #ggtitle('Percentage of valid observations of KdPAR (2013-2022)') +
  ylab('KdPAR (/m)') + xlab('Radius (km)') +
  theme_light() +
  theme(legend.position = 'bottom')
figS2

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/FigureS2_KdPAR_boxplot_CLASS2_vs_Radii.png',width = 8, height = 6, dpi=600, figS2)

3.3 Fig S3 KdPAR 90th percentile vs Catchment Size

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile90' & CATCHMT_AR>0 & CLASS2 == 'Open Coast') %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  ylab('KdPAR 10yr (2013-2022) 90th percentile (/m)') + xlab('Catchment Size (km2)') +
  ggtitle('Open Coast') +
  theme_light() +
  theme(legend.position = 'bottom')
p

`geom_smooth()` using formula = 'y ~ x'

p2 <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & Metrics == 'percentile90' & CATCHMT_AR>0 & CLASS2 == 'Enclosed Coast') %>% 
  ggplot(data=.,aes(x=CATCHMT_AR,y=KdPAR)) +
  geom_point() +
  geom_smooth(method = 'lm') +
  ylab('KdPAR 10yr (2013-2022) 90th percentile (/m)') + xlab('Catchment Size (km2)') +
  ggtitle('Enclosed Coast') +
  theme_light() +
  theme(legend.position = 'bottom')
p2

`geom_smooth()` using formula = 'y ~ x'

figS3 <- ggarrange(p,
                     p2,
                     labels = c("A", "B"),
                     ncol = 2, nrow = 1)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

figS3

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/FigureS3_KdPAR_90th_CatchmentSize.png',width = 8, height = 6, dpi=600, figS3)

3.4 Fig S4 Same as Figure 8 but for coastal council regions

p <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & CLASS2 == "Open Coast" & Metrics %in% c('percentile10','median','percentile90')) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
    mutate(Metrics = factor(Metrics,levels = c('percentile5','10th percentile','Median','90th percentile','percentile95'))) %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~Metrics,ncol=3) +
  scale_fill_viridis_d() +
  ylab('KdPAR (/m)') + xlab('Bioregions') +
  ggtitle("Open Coast") +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p

p2 <- kd_metrics %>% 
  mutate(Buffer_Size = Buffer_Size/1000) %>% 
  filter(Buffer_Size == 10 & !is.na(KdPAR) & CLASS2 == "Enclosed Coast" & Metrics %in% c('percentile10','median','percentile90')) %>% 
  mutate(Metrics = case_when(Metrics == 'percentile10' ~  '10th percentile',
                             Metrics == 'median' ~  'Median',
                             Metrics == 'percentile90' ~  '90th percentile')) %>%
    mutate(Metrics = factor(Metrics,levels = c('percentile5','10th percentile','Median','90th percentile','percentile95'))) %>% 
  ggplot(data=.,aes(x=factor(REGC2017_N),y=KdPAR)) +
  geom_boxplot() +
  facet_wrap(~Metrics,ncol=3) +
  scale_fill_viridis_d() +
  ylab('KdPAR (/m)') + xlab('Bioregions') +
  ggtitle("Enclosed Coast") +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
p2

figureS4 <- ggarrange(p,
                     p2,
                     labels = c("A", "B"),
                     ncol = 1, nrow = 2)
figureS4

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/FigureS4_Kd_vs_CouncilRegions.png',width = 12, height = 8, dpi=600, figureS4)

3.5 Figure S5 Number of valid observations per class2, 10 km radius

p <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'valid_days_pc') %>% 
  ggplot(data=.) +
  geom_boxplot(aes(x=factor(CLASS2),y=KdPAR)) +
  scale_fill_viridis_d() +
  ggtitle("500 m KdPAR product") +
  ylab('Percent of valid images (%)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

p2 <- kd_metrics %>% 
  filter(Buffer_Size == 10000 & Metrics == 'valid_days_pc') %>% 
  ggplot(data=.) +
  geom_boxplot(aes(x=factor(CLASS2),y=Kd490)) +
  scale_fill_viridis_d() +
  ggtitle("4 km Kd490 product") +
  ylab('Percent of valid images (%)') + xlab('') +
  theme_light() +
  theme(legend.position = 'bottom',
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

figureS5 <- ggarrange(p2,
                     p,
                     labels = c("A", "B"),
                     ncol = 2, nrow = 1)
figureS5

#ggsave('C:/Users/thoralf/OneDrive - NIWA/Documents/Postdoc_TauKiAakau/Goal_3_LandUse_Effects_WaterClarity/Figures/Feb2025/FigureS5_PercentValidObs_Kd490_KdPAR_10km.png',width = 8, height = 6, dpi=600, figureS5)

References

Gall, Mark P, Matthew H Pinkerton, Tilmann Steinmetz, and Simon Wood. 2022. “Satellite Remote Sensing of Coastal Water Quality in New Zealand.” New Zealand Journal of Marine and Freshwater Research 56 (3): 585–616.