psi_model_analysis

library(dplyr)
library(tidyverse)
library(lubridate)
library(stringr)
library(kableExtra)
library(ggplot2)
library(gridExtra)
library(multcompView)

dat = read.csv("C://Users/bdmor/Box/NASA-IAV/classification/s2_psi_analysis/outputs/psi_model_stats/psi_model_outputs.csv")

dat$date = as_date(dat$date, tz = NULL)
dat$band_combo = paste(dat$band_y, dat$band_x, sep = "_")

How many models were Forward or Reverse?

kable(as.data.frame(table(dat$version))) %>%
  kable_material()

Var1	Freq
EITHER	42
FORWARD	189
REVERSE	3381

Reverse is the most common version of PSI calculation Get rid of models that ran FORWARD

x = which(dat$version == "FORWARD")
dat = dat[-x,]

Only keep bands combo’s that have the same resolution (e.g., no mixing 10m x 20m)

x = which(dat$band_x == "B8A")
dat = dat[x,]

dat = dat[complete.cases(dat),]

# band_combos = expand.grid(c('B2', 'B3', 'B4'), c('B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12'), stringsAsFactors = FALSE)
# band_combos = paste(band_combos$Var1, band_combos$Var2, sep = "_")

dat$date = as_date(dat$date)
dat$month = month(dat$date)#, label = TRUE, abbr = TRUE)
dat$year = year(dat$date)

Make seasons category

dat$season = NA
x = which(dat$month >2 & dat$month < 6)
dat$season[x] = "MAM"
x = which(dat$month >5 & dat$month < 9)
dat$season[x] = "JJA"
x = which(dat$month >8 & dat$month < 12)
dat$season[x] = "SON"
x = which(dat$month < 3 | dat$month == 12)
dat$season[x] = "DJF"

Summary table of model run statistics

df = dat %>%
  select(band_combo, del_acc)%>%
  group_by(band_combo) %>%
  summarise(
    min = min(del_acc),
    Q25 = quantile(del_acc, 0.25),
    Q50 = quantile(del_acc, 0.5),
    mean = mean(del_acc),
    Q75 = quantile(del_acc, 0.75),
    max = max(del_acc),
    sd = sd(del_acc)
  ) %>%
  arrange(desc(Q50))

track = data.frame(rank = 1:3)
track$del_acc = df$band_combo

sum.df = data.frame(band_combo = unique(dat$band_combo))
sum.df$del_acc[match(sum.df$band_combo, df$band_combo)] = df$Q50

Delta Accuracy Summary

df = dat %>%
  select(band_combo, del_kap)%>%
  group_by(band_combo) %>%
  summarise(
    min = min(del_kap),
    Q25 = quantile(del_kap, 0.25),
    Q50 = quantile(del_kap, 0.5),
    mean = mean(del_kap),
    Q75 = quantile(del_kap, 0.75),
    max = max(del_kap),
    sd = sd(del_kap)
  ) %>%
  arrange(desc(Q50))

track$del_kap = df$band_combo

sum.df$del_kap[match(sum.df$band_combo, df$band_combo)] = df$Q50

Delta Kappa Summary

df = dat %>%
  select(band_combo, sav_acc)%>%
  group_by(band_combo) %>%
  summarise(
    min = min(sav_acc),
    Q25 = quantile(sav_acc, 0.25),
    Q50 = quantile(sav_acc, 0.5),
    mean = mean(sav_acc),
    Q75 = quantile(sav_acc, 0.75),
    max = max(sav_acc),
    sd = sd(sav_acc)
  ) %>%
  arrange(desc(Q50))

track$sav_acc = df$band_combo

sum.df$sav_acc[match(sum.df$band_combo, df$band_combo)] = df$Q50

SAV Accuracy Summary

df = dat %>%
  select(band_combo, sav_kap)%>%
  group_by(band_combo) %>%
  summarise(
    min = min(sav_kap),
    Q25 = quantile(sav_kap, 0.25),
    Q50 = quantile(sav_kap, 0.5),
    mean = mean(sav_kap),
    Q75 = quantile(sav_kap, 0.75),
    max = max(sav_kap),
    sd = sd(sav_kap)
  ) %>%
  arrange(desc(Q50))

track$sav_kap = df$band_combo

sum.df$sav_kap[match(sum.df$band_combo, df$band_combo)] = df$Q50

SAV Kappa Summary

df = dat %>%
  select(band_combo, wat_acc)%>%
  group_by(band_combo) %>%
  summarise(
    min = min(wat_acc),
    Q25 = quantile(wat_acc, 0.25),
    Q50 = quantile(wat_acc, 0.5),
    mean = mean(wat_acc),
    Q75 = quantile(wat_acc, 0.75),
    max = max(wat_acc),
    sd = sd(wat_acc)
  ) %>%
  arrange(desc(Q50))

track$wat_acc = df$band_combo

sum.df$wat_acc[match(sum.df$band_combo, df$band_combo)] = df$Q50

Water Accuracy Summary

df = dat %>%
  select(band_combo, wat_kap)%>%
  group_by(band_combo) %>%
  summarise(
    min = min(wat_kap),
    Q25 = quantile(wat_kap, 0.25),
    Q50 = quantile(wat_kap, 0.5),
    mean = mean(wat_kap),
    Q75 = quantile(wat_kap, 0.75),
    max = max(wat_kap),
    sd = sd(wat_kap)
  ) %>%
  arrange(desc(Q50))

track$wat_kap = df$band_combo

sum.df$wat_kap[match(sum.df$band_combo, df$band_combo)] = df$Q50

Water Kappa Summary

Band Combination with B4 and B8A had highest PSI Delta Accuracy, Kappa, SAV Accuracy and Kappa. Only B2/B8A had higher Water Accuracy and Kappa.

kable(sum.df) %>%
  kable_material()

band_combo	del_acc	del_kap	sav_acc	sav_kap	wat_acc	wat_kap
B4_B8A	70.000	0.4000	64.2575	0.285	50.8330	0.017
B3_B8A	64.705	0.2860	63.0250	0.261	52.5000	0.050
B2_B8A	68.200	0.3665	64.1670	0.283	51.3945	0.028

rank.vect = unlist(track[1, 2:ncol(track)])
rank.df = as.data.frame(table(rank.vect))
names(rank.df) = c("band_combo", "Freq")

kable(rank.df) %>%
  kable_material()

band_combo	Freq
B2_B8A	2
B4_B8A	4

Summary of PSI Accuracy and Kappa

df = dat %>%
  select(band_combo, del_acc)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(del_acc ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

max.acc = max(df$del_acc)

psi.acc.plot = ggplot(df, aes(x = band_combo, y = del_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = max.acc+5, label = label), size = 5) +
  labs(x = "Band Combination", y = "PSI Accuracy", title = "Box Plot of PSI Delta Accuracy by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, del_kap)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(del_kap ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

max.kap = max(df$del_kap)

psi.kap.plot = ggplot(df, aes(x = band_combo, y = del_kap))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = max.kap+0.05, label = label), size = 5) +
  labs(x = "Band Combination", y = "PSI Kappa", title = "Box Plot of PSI Delta Kappa by Band Combination")+
  theme_classic()

Box plots of PSI Delta Accuracy and Kappa

grid.arrange(psi.acc.plot, psi.kap.plot, nrow = 2)

Summary of PSI SAV Model Stats

df = dat %>%
  select(band_combo, sav_pval)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(sav_pval ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)


sav.p.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, sav_pval) %>%
  summarise(
    min = boxplot.stats(sav_pval)$stats[1],
    max = boxplot.stats(sav_pval)$stats[5]
  )

df.nol = df %>%
  inner_join(sav.p.ol, by = 'band_combo')%>%
  filter(sav_pval > min, sav_pval < max) %>%
  select(band_combo, sav_pval)

max.val = max(sav.p.ol$max)
  

sav.p.plot =  ggplot(df.nol, aes(x = band_combo, y = sav_pval))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 1.05, label = label), size = 5) +
  labs(x = "Band Combination", y = "p-Value", title = "Box Plot of PSI SAV p-Value by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, sav_rsq)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(sav_rsq ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$sav_rsq)
sav.rsq.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, sav_rsq) %>%
  summarise(
    min = boxplot.stats(sav_rsq)$stats[1],
    max = boxplot.stats(sav_rsq)$stats[5]
  )

df.nol = df %>%
  inner_join(sav.rsq.ol, by = 'band_combo')%>%
  filter(sav_rsq > min, sav_rsq < max) %>%
  select(band_combo, sav_rsq)

max.val = max(sav.rsq.ol$max)

sav.r2.plot = ggplot(df.nol, aes(x = band_combo, y = sav_rsq))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 1, label = label), size = 5) +
  labs(x = "Band Combination", y = "R2", title = "Box Plot of PSI SAV R2 by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, sav_slp)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(sav_slp ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$sav_slp)
sav.slp.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, sav_slp) %>%
  summarise(
    min = boxplot.stats(sav_slp)$stats[1],
    max = boxplot.stats(sav_slp)$stats[5]
  )

df.nol = df %>%
  inner_join(sav.slp.ol, by = 'band_combo')%>%
  filter(sav_slp > min, sav_slp < max) %>%
  select(band_combo, sav_slp)

max.val = round(max(sav.slp.ol$max), digits = 2)

sav.slp.plot = ggplot(df.nol, aes(x = band_combo, y = sav_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = max.val, label = label), size = 5) +
  labs(x = "Band Combination", y = "Slope", title = "Box Plot of PSI SAV Slope by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, sav_int)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(sav_int ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)


sav.int.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, sav_int) %>%
  summarise(
    min = boxplot.stats(sav_int)$stats[1],
    max = boxplot.stats(sav_int)$stats[5]
  )

df.nol = df %>%
  inner_join(sav.int.ol, by = 'band_combo')%>%
  filter(sav_int > min, sav_int < max) %>%
  select(band_combo, sav_int)

max.val = round(max(sav.int.ol$max), digits = 3)
  

sav.int.plot =  ggplot(df.nol, aes(x = band_combo, y = sav_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.065, label = label), size = 5) +
  labs(x = "Band Combination", y = "Intercept", title = "Box Plot of PSI SAV Intercept by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, sav_acc)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(sav_acc ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$sav_acc)
sav.acc.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, sav_acc) %>%
  summarise(
    min = boxplot.stats(sav_acc)$stats[1],
    max = boxplot.stats(sav_acc)$stats[5]
  )

df.nol = df %>%
  inner_join(sav.acc.ol, by = 'band_combo')%>%
  filter(sav_acc > min, sav_acc < max) %>%
  select(band_combo, sav_acc)

max.val = max(sav.acc.ol$max)


sav.acc.plot = ggplot(df.nol, aes(x = band_combo, y = sav_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Band Combination", y = "Accuracy", title = "Box Plot of PSI SAV Accuracy by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, sav_kap)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(sav_kap ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$sav_kap)
sav.kap.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, sav_kap) %>%
  summarise(
    min = boxplot.stats(sav_kap)$stats[1],
    max = boxplot.stats(sav_kap)$stats[5]
  )

df.nol = df %>%
  inner_join(sav.kap.ol, by = 'band_combo')%>%
  filter(sav_kap > min, sav_kap < max) %>%
  select(band_combo, sav_kap)

max.val = max(sav.kap.ol$max)

sav.kap.plot = ggplot(df.nol, aes(x = band_combo, y = sav_kap))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = max.val, label = label), size = 5) +
  labs(x = "Band Combination", y = "Kappa", title = "Box Plot of PSI SAV Kappa by Band Combination")+
  theme_classic()

Box plots of PSI SAV Models

Note: Some plots have a wide range of outliers that make viewing the boxplots difficult. The most Extreme outliers have been removed for plotting purposes, but Tukey test is still on all data with outliers.

grid.arrange(sav.p.plot, sav.r2.plot, sav.slp.plot, sav.int.plot, sav.acc.plot, sav.kap.plot, nrow = 6)

Summary of PSI Water Model Stats

df = dat %>%
  select(band_combo, wat_pval)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(wat_pval ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$wat_pval)
wat.p.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, wat_pval) %>%
  summarise(
    min = boxplot.stats(wat_pval)$stats[1],
    max = boxplot.stats(wat_pval)$stats[5]
  )

df.nol = df %>%
  inner_join(wat.p.ol, by = 'band_combo')%>%
  filter(wat_pval > min, wat_pval < max) %>%
  select(band_combo, wat_pval)

max.val = max(wat.p.ol$max)

wat.p.plot = ggplot(df.nol, aes(x = band_combo, y = wat_pval))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.015, label = label), size = 5) +
  labs(x = "Band Combination", y = "p-Value", title = "Box Plot of PSI Water p-Value by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, wat_rsq)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(wat_rsq ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$wat_rsq)
wat.rsq.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, wat_rsq) %>%
  summarise(
    min = boxplot.stats(wat_rsq)$stats[1],
    max = boxplot.stats(wat_rsq)$stats[5]
  )

df.nol = df %>%
  inner_join(wat.rsq.ol, by = 'band_combo')%>%
  filter(wat_rsq > min, wat_rsq < max) %>%
  select(band_combo, wat_rsq)

max.val = max(wat.rsq.ol$max)

wat.r2.plot = ggplot(df.nol, aes(x = band_combo, y = wat_rsq))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 1, label = label), size = 5) +
  labs(x = "Band Combination", y = "R2", title = "Box Plot of PSI Water R2 by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, wat_slp)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(wat_slp ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$wat_slp)
wat.slp.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, wat_slp) %>%
  summarise(
    min = boxplot.stats(wat_slp)$stats[1],
    max = boxplot.stats(wat_slp)$stats[5]
  )

df.nol = df %>%
  inner_join(wat.slp.ol, by = 'band_combo')%>%
  filter(wat_slp > min, wat_slp < max) %>%
  select(band_combo, wat_slp)

max.val = max(wat.slp.ol$max)

wat.slp.plot = ggplot(df.nol, aes(x = band_combo, y = wat_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 2.38, label = label), size = 5) +
  labs(x = "Band Combination", y = "Slope", title = "Box Plot of PSI Water Slope by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, wat_int)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(wat_int ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$wat_int)
wat.int.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, wat_int) %>%
  summarise(
    min = boxplot.stats(wat_int)$stats[1],
    max = boxplot.stats(wat_int)$stats[5]
  )

df.nol = df %>%
  inner_join(wat.int.ol, by = 'band_combo')%>%
  filter(wat_int > min, wat_int < max) %>%
  select(band_combo, wat_int)

max.val = max(wat.int.ol$max)

wat.int.plot = ggplot(df.nol, aes(x = band_combo, y = wat_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.06, label = label), size = 5) +
  labs(x = "Band Combination", y = "Intercept", title = "Box Plot of PSI Water Intercept by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, wat_acc)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(wat_acc ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$wat_acc)
wat.acc.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, wat_acc) %>%
  summarise(
    min = boxplot.stats(wat_acc)$stats[1],
    max = boxplot.stats(wat_acc)$stats[5]
  )

df.nol = df %>%
  inner_join(wat.acc.ol, by = 'band_combo')%>%
  filter(wat_acc > min, wat_acc < max) %>%
  select(band_combo, wat_acc)

max.val = max(wat.acc.ol$max)

wat.acc.plot = ggplot(df.nol, aes(x = band_combo, y = wat_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 63, label = label), size = 5) +
  labs(x = "Band Combination", y = "Accuracy", title = "Box Plot of PSI Water Accuracy by Band Combination")+
  theme_classic()

df = dat %>%
  select(band_combo, wat_kap)

comps = combn(unique(df$band_combo), 2, simplify = TRUE) 
comps = asplit(comps, 2)

amodel = aov(wat_kap ~ band_combo, df)
tukey_result <- TukeyHSD(amodel, ordered = T)
ps = tukey_result$band_combo[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

# max.acc = max(df$wat_kap)
wat.kap.ol = df %>%
  group_by(band_combo)%>%
  select(band_combo, wat_kap) %>%
  summarise(
    min = boxplot.stats(wat_kap)$stats[1],
    max = boxplot.stats(wat_kap)$stats[5]
  )

df.nol = df %>%
  inner_join(wat.kap.ol, by = 'band_combo')%>%
  filter(wat_kap > min, wat_kap < max) %>%
  select(band_combo, wat_kap)

max.val = max(wat.kap.ol$max)

wat.kap.plot = ggplot(df.nol, aes(x = band_combo, y = wat_kap))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = max.val, label = label), size = 5) +
  labs(x = "Band Combination", y = "Kappa", title = "Box Plot of PSI Water Kappa by Band Combination")+
  theme_classic()

Box plots of PSI Water Models

Note: Some plots have a wide range of outliers that make viewing the boxplots difficult. The most Extreme outliers have been removed for plotting purposes, but Tukey test is still on all data with outliers.

grid.arrange(wat.p.plot, wat.r2.plot, wat.slp.plot, wat.int.plot, wat.acc.plot, wat.kap.plot, nrow = 6)

Histogram plots of PSI SAV vs. PSI Water Slopes

df = dat %>%
  select(band_combo, sav_slp, wat_slp) %>%
  pivot_longer(cols = contains("slp"),
               names_to = "class",
               values_to = "slp") %>%
    mutate(class = str_replace(
      class, "_slp", ""))

ol = boxplot.stats(df$slp)

df.nol = df %>%
  filter(slp > ol$stats[1], slp < ol$stats[5]) %>%
  select(band_combo, class, slp)

    
ggplot(df.nol, aes(x = slp, group = class, fill = class))+
  geom_histogram(position = "identity", binwidth = 0.1, color = 'black', alpha = 0.5)+
  scale_fill_manual(labels = c("sav", "wat"), values = c("chartreuse4", "blue"))+
  theme_classic()+
  labs(x = "Slope", y = "Count", title = "Histograms of SAV vs. Wat PSI Slopes by Band Combinations", group = "Class", fill = "Class")+
  facet_wrap(~band_combo, nrow = 3)

Does PSI differ over time?

df = dat %>%
  filter(band_combo == "B4_B8A") 
df$month = month(df$date, label = TRUE, abbr = TRUE)
df$year = as.factor(df$year)
df$season = factor(df$season, labels = c("DJF", "MAM", "JJA", "SON"), ordered = TRUE)

# comps_year = combn(unique(year(df$date)), 2, simplify = TRUE)
# comps_year = asplit(comps_year, 2)
# 
# comps_month = combn(unique(month(df$date)), 2, simplify = TRUE)
# comps_month = asplit(comps_month, 2)
# 
# comps_season = combn(unique(df$season), 2, simplify = TRUE)
# comps_season = asplit(comps_season, 2)

Boxplots of Delta PSI Accuracy

# Del accuracy
amodel = aov(del_acc ~ year, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$year[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

del_acc_p1 = ggplot(df, aes(x = year, y = del_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Year", y = "Accuracy", title = "Year")+
  theme_classic()


amodel = aov(del_acc ~ season, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$season[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

del_acc_p2 = ggplot(df, aes(x = season, y = del_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Season", y = "Accuracy", title = "Season")+
  theme_classic()


amodel = aov(del_acc ~ month, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$month[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

del_acc_p3 = ggplot(df, aes(x = month, y = del_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Month", y = "Accuracy", title = "Month")+
  theme_classic()

grid.arrange(del_acc_p1, del_acc_p2, del_acc_p3, nrow = 3,
             top = "PSI Delta Accuracy Temporal Trends B4-B8A")

Boxplots of SAV Accuracy

# sav acc
amodel = aov(sav_acc ~ year, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$year[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_acc_p1 = ggplot(df, aes(x = year, y = sav_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Year", y = "Accuracy", title = "Year")+
  theme_classic()


amodel = aov(sav_acc ~ season, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$season[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_acc_p2 = ggplot(df, aes(x = season, y = sav_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Season", y = "Accuracy", title = "Season")+
  theme_classic()


amodel = aov(sav_acc ~ month, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$month[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_acc_p3 = ggplot(df, aes(x = month, y = sav_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Month", y = "Accuracy", title = "Month")+
  theme_classic()

grid.arrange(sav_acc_p1, sav_acc_p2, sav_acc_p3, nrow = 3,
             top = "SAV Accuracy Temporal Trends B4-B8A")

Boxplots of Water Accuracy

# wat acc
amodel = aov(wat_acc ~ year, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$year[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_acc_p1 = ggplot(df, aes(x = year, y = wat_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Year", y = "Accuracy", title = "Year")+
  theme_classic()


amodel = aov(wat_acc ~ season, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$season[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_acc_p2 = ggplot(df, aes(x = season, y = wat_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Season", y = "Accuracy", title = "Season")+
  theme_classic()


amodel = aov(wat_acc ~ month, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$month[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_acc_p3 = ggplot(df, aes(x = month, y = wat_acc))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 100, label = label), size = 5) +
  labs(x = "Month", y = "Accuracy", title = "Month")+
  theme_classic()

grid.arrange(wat_acc_p1, wat_acc_p2, wat_acc_p3, nrow = 3,
             top = "Water Accuracy Temporal Trends B4-B8A")

Boxplots of SAV Slope

# sav slp
amodel = aov(sav_slp ~ year, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$year[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_slp_p1 = ggplot(df, aes(x = year, y = sav_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 2.85, label = label), size = 5) +
  labs(x = "Year", y = "Slope", title = "Year")+
  theme_classic()


amodel = aov(sav_slp ~ season, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$season[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_slp_p2 = ggplot(df, aes(x = season, y = sav_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 2.85, label = label), size = 5) +
  labs(x = "Season", y = "Slope", title = "Season")+
  theme_classic()


amodel = aov(sav_slp ~ month, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$month[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_slp_p3 = ggplot(df, aes(x = month, y = sav_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 2.85, label = label), size = 5) +
  labs(x = "Month", y = "Slope", title = "Month")+
  theme_classic()

grid.arrange(sav_slp_p1, sav_slp_p2, sav_slp_p3, nrow = 3,
             top = "SAV Slope Temporal Trends B4-B8A")

Boxplots of Water Slope

# wat slp
amodel = aov(wat_slp ~ year, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$year[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_slp_p1 = ggplot(df, aes(x = year, y = wat_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 4.82, label = label), size = 5) +
  labs(x = "Year", y = "Slope", title = "Year")+
  theme_classic()


amodel = aov(wat_slp ~ season, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$season[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_slp_p2 = ggplot(df, aes(x = season, y = wat_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 4.82, label = label), size = 5) +
  labs(x = "Season", y = "Slope", title = "Season")+
  theme_classic()


amodel = aov(wat_slp ~ month, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$month[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_slp_p3 = ggplot(df, aes(x = month, y = wat_slp))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 4.82, label = label), size = 5) +
  labs(x = "Month", y = "Slope", title = "Month")+
  theme_classic()

grid.arrange(wat_slp_p1, wat_slp_p2, wat_slp_p3, nrow = 3,
             top = "Water Slope Temporal Trends B4-B8A")

Boxplots of SAV Intercept

# sav int
amodel = aov(sav_int ~ year, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$year[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_int_p1 = ggplot(df, aes(x = year, y = sav_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.072, label = label), size = 5) +
  labs(x = "Year", y = "Intercept", title = "Year")+
  theme_classic()


amodel = aov(sav_int ~ season, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$season[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_int_p2 = ggplot(df, aes(x = season, y = sav_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.072, label = label), size = 5) +
  labs(x = "Season", y = "Intercept", title = "Season")+
  theme_classic()


amodel = aov(sav_int ~ month, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$month[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

sav_int_p3 = ggplot(df, aes(x = month, y = sav_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.072, label = label), size = 5) +
  labs(x = "Month", y = "Intercept", title = "Month")+
  theme_classic()

grid.arrange(sav_int_p1, sav_int_p2, sav_int_p3, nrow = 3,
             top = "SAV Intercept Temporal Trends B4-B8A")

Boxplots of Water Intercept

# wat int
amodel = aov(wat_int ~ year, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$year[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_int_p1 = ggplot(df, aes(x = year, y = wat_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.071, label = label), size = 5) +
  labs(x = "Year", y = "Intercept", title = "Year")+
  theme_classic()


amodel = aov(wat_int ~ season, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$season[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_int_p2 = ggplot(df, aes(x = season, y = wat_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.071, label = label), size = 5) +
  labs(x = "Season", y = "Intercept", title = "Season")+
  theme_classic()


amodel = aov(wat_int ~ month, df)
tukey_result <- TukeyHSD(amodel)
ps = tukey_result$month[, 'p adj']
cld <- multcompLetters(ps)

cld_data <- data.frame(
  group = names(cld$Letters),
  label = cld$Letters
)

wat_int_p3 = ggplot(df, aes(x = month, y = wat_int))+
  geom_boxplot()+
  geom_text(data = cld_data, aes(x = group, y = 0.071, label = label), size = 5) +
  labs(x = "Month", y = "Intercept", title = "Month")+
  theme_classic()

grid.arrange(wat_int_p1, wat_int_p2, wat_int_p3, nrow = 3,
             top = "Water Intercept Temporal Trends B4-B8A")