library("tidyr")
library("dplyr")
library("readr")
library("ggplot2")
library("ggbreak")
library("patchwork")
library("ggpubr")

MyMerge       <- function(x, y){
  df            <- merge(x, y, all = TRUE)
  rownames(df)  <- df$Row.names
  df$Row.names  <- NULL
  return(df)
}

Team vole NDVI/COVER

Load in NDVI data

TV_NDVI_orig_2021 <- read.csv("../Data/2021_Toolik_NDVI.csv", header = TRUE) 
TV_NDVI_orig_2020 <- read.csv("../Data/2020_Toolik_NDVI.csv", header = TRUE) 
TV_NDVI_orig_2019 <- read.csv("../Data/2019_Toolik_NDVI.csv", header = TRUE) 
TV_NDVI_orig_2018 <- read.csv("../Data/2018_Toolik_NDVI.csv", header = TRUE)

Prep the Team Vole NDVI data

Remove data we don’t need

TV_NDVI_clean_2021 <-  TV_NDVI_orig_2021 %>%
  select(Site, Treatment, Quadrat.Chamber.., Replicate.., NDVI) %>%
  filter(Quadrat.Chamber.. != "CH1" & Quadrat.Chamber.. != "CH2" & Quadrat.Chamber.. != "CH3" & Quadrat.Chamber.. != "") %>%
  rename("quad" = "Quadrat.Chamber..", "treatment" = "Treatment", "site" = "Site", "rep" = "Replicate..") %>%
  mutate(year = "2021")
    
      TV_NDVI_clean_2021$quad <- TV_NDVI_clean_2021$quad %>%
        recode("Q1" = "1", "Q2"= "2", "Q3"= "3", "Q4"= "4", "Q5"= "5", "Q6"= "6", "Q7"= "7", "Q8"= "8") 


TV_NDVI_clean_2020 <-  TV_NDVI_orig_2020 %>%
  select(Site, Treatment, Quadrat, Scan, NDVI) %>%
  rename("quad" = "Quadrat", "treatment" = "Treatment", "site" = "Site", "rep" = "Scan")%>%
  filter(quad != "MISCAN")%>%
  mutate(year = "2020")

TV_NDVI_clean_2019 <-  TV_NDVI_orig_2019 %>%
  select(Site, Treatment, Quadrat, NDVI) %>%
  filter(Quadrat != "CH1" & Quadrat != "CH2" & Quadrat != "CH3") %>%
  rename("quad" = "Quadrat", "treatment" = "Treatment", "site" = "Site")%>%
  mutate(year = "2019")

      TV_NDVI_clean_2019$quad <- TV_NDVI_clean_2019$quad %>%
        recode("Q1" = "1", "Q2"= "2", "Q3"= "3", "Q4"= "4", "Q5"= "5", "Q6"= "6", "Q7"= "7", "Q8"= "8")

TV_NDVI_clean_2018 <-  TV_NDVI_orig_2018 %>%
  select(SITE, TREATMENT, QUAD2, NDVI) %>%
  filter(QUAD2 != "F1" & QUAD2 != "F2" & QUAD2 != "F3" & QUAD2 != "") %>%
  rename("quad" = "QUAD2", "treatment" = "TREATMENT", "site" = "SITE")%>%
  mutate(year = "2018")

      TV_NDVI_clean_2018$quad <- TV_NDVI_clean_2018$quad %>%
       recode("C1" = "1", "C2"= "2", "C3"= "3", "C4"= "4", "C5"= "5", "C6"= "6", "C7"= "7", "C8"= "8")
      
unique(TV_NDVI_clean_2018$quad)

## [1] "1" "2" "3" "4" "5" "6" "7" "8"

Merge the years together

TV_NDVI_CLEAN           <- Reduce(MyMerge, list(TV_NDVI_clean_2021, TV_NDVI_clean_2020, TV_NDVI_clean_2019, TV_NDVI_clean_2018))%>%
  mutate(quad = as.numeric(quad))

Average across reps, then across quadrants.

NOTE: 2018 and 2019 didn’t have columns to indicate rep number, they just listed the same quad info multiple times. This will still work with the code below because we do not need to group by rep.
LOGIC CHECK: N should equal 8 since there should have been 8 reps. 7 is also acceptable (we might have only done 7 or we could have removed a scan for an error). If any other numbers show up then investigate why its not 7 or 8 and correct the mistake.

TV_NDVI_avg_rep <- TV_NDVI_CLEAN %>%
  group_by(year, site, treatment, quad) %>%
  summarise(NDVI_avg = mean(NDVI, na.rm = TRUE), n = n(), sd_rep = sd(NDVI, na.rm = TRUE))

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

TV_NDVI_avg <- TV_NDVI_avg_rep %>%
  group_by(year, site, treatment) %>%
  summarise(avg_NDVI = mean(NDVI_avg, na.rm = TRUE), n = n(), sd = sd(NDVI_avg, na.rm = TRUE))

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

Box plot - NDVI

TV_NDVI_avg$treatment <- factor(TV_NDVI_avg$treatment, levels = c("CT", "EX", "PU", "PR"))

  ggplot(TV_NDVI_avg, aes(x= year, y=avg_NDVI, fill = treatment))+
    geom_hline(yintercept = 0.6, color="#D0CECE")+
      geom_boxplot(position=position_dodge(.85))+
    scale_fill_manual(values=c("CT" = "#a6611a", 
                                "EX" = "#dfc27d",
                                "PU" = "#80cdc1",
                                "PR" = "#018571"))+
  stat_summary(fun.y="mean", color="#D1D1D1", size = .5, shape = 18, position=position_dodge(.85))+
     theme_light()+
     theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
     theme(axis.title.y = element_text(face = "bold"))+
     theme(axis.title.x = element_text(face = "bold"))+
     theme(plot.title = element_text(hjust = .5))+
     theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
     theme(axis.text.y = element_text(color = "black", size = 8))+
     theme(aspect.ratio = 9/18.5)+
  labs(y = "Mean NDVI", x = "")

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 15 rows containing missing values (geom_segment).

    #theme(legend.position="bottom")
  
ggsave("../Figures/Teamvole_boxplot_NDVI.jpeg")

## Saving 7 x 5 in image

## Warning: Removed 15 rows containing missing values (geom_segment).

TV_NDVI_avg_rep$treatment <- factor(TV_NDVI_avg_rep$treatment, levels = c("CT", "EX", "PU", "PR"))

  ggplot(TV_NDVI_avg_rep, aes(x= year, y=NDVI_avg, fill = treatment))+
    geom_hline(yintercept = 0.6, color="#D0CECE")+
      geom_boxplot(position=position_dodge(.85))+
    scale_fill_manual(values=c("CT" = "#a6611a", 
                                "EX" = "#dfc27d",
                                "PU" = "#80cdc1",
                                "PR" = "#018571"))+
    stat_summary(fun.y="mean", color="#D1D1D1", size = .1, shape = 18, position=position_dodge(.85))+
     theme_light()+
     theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
     theme(axis.title.y = element_text(face = "bold"))+
     theme(axis.title.x = element_text(face = "bold"))+
     theme(plot.title = element_text(hjust = .5))+
     theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
     theme(axis.text.y = element_text(color = "black", size = 8))+
     theme(aspect.ratio = 9/18.5)+
  labs(y = "Mean NDVI", x = "")+
    #theme(legend.position="bottom")+
     facet_grid(site ~ .)

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 15 rows containing missing values (geom_segment).

## Warning: Removed 15 rows containing missing values (geom_segment).

## Warning: Removed 15 rows containing missing values (geom_segment).

ggsave("../Figures/Teamvole_boxplot_NDVI_sites.jpeg")

## Saving 7 x 5 in image

## Warning: Removed 15 rows containing missing values (geom_segment).

## Warning: Removed 15 rows containing missing values (geom_segment).

## Warning: Removed 15 rows containing missing values (geom_segment).

Load in Team Vole cover data

TV_cover_orig <- read.csv("../Derived CSVs/TeamVole_relcov_funcgrp_ALL.csv", header = TRUE)

Prep TV cover

Only looking at a few species, so the rest can be deleted
Need to change cover from long to wide format so it can be plot against NDVI data

TV_cover_subset <- TV_cover_orig %>%
  select(year, site, treatment, quad, species, cover)%>%
  filter(species == "Eri vag" | species == "Bet nan" | species == "Rub cha") %>%
  pivot_wider(names_from = species, values_from = cover) %>%
  rename("Bet_nan" = "Bet nan", "Eri_vag" = "Eri vag", "Rub_cha" = "Rub cha") %>%
  mutate(year = as.character(year), quad = as.numeric(quad))

Join NDVI and COVER data

Join the dataset, drop columns we don’t need, delete 2020 data as we don’t have cover for that year

TV_join <- left_join(TV_NDVI_avg_rep, TV_cover_subset, by= c("year", "site", "treatment", "quad")) %>%
  select(-n, sd_rep) %>%
  filter(year != "2020")

Scatter plot of NDVI vs species cover

Bet nan vs NDVI

highlight_CT <- TV_join %>% filter(treatment == "CT")
highlight_EX <- TV_join %>% filter(treatment == "EX")
highlight_PU <- TV_join %>% filter(treatment == "PU")
highlight_PR <- TV_join %>% filter(treatment == "PR")

ggplot(data = TV_join, aes(x = Bet_nan, y = NDVI_avg)) +
  geom_point(size = 2)+
  geom_smooth(method = "lm", se = TRUE, color = "black")+
   geom_point(data=highlight_CT, aes(x=Bet_nan,y=NDVI_avg), size = 2, color = "#a6611a")+
  geom_point(data=highlight_EX, aes(x=Bet_nan,y=NDVI_avg), size = 2, color = "#dfc27d")+
  geom_point(data=highlight_PU, aes(x=Bet_nan,y=NDVI_avg), size = 2, color = "#80cdc1")+
  geom_point(data=highlight_PR, aes(x=Bet_nan,y=NDVI_avg), size = 2, color = "#018571")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Bet nan")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.3)+ #adds R2 
  #stat_cor(method = "pearson", label.x = .3)+ #adds p-value 
  facet_grid(site ~ .)

  #facet_grid(. ~ site)

#ggsave("../Figures/Teamvole_regrestion_Betnan_NDVI.jpeg")

Eri vag vs NDVI

highlight_CT <- TV_join %>% filter(treatment == "CT")
highlight_EX <- TV_join %>% filter(treatment == "EX")
highlight_PU <- TV_join %>% filter(treatment == "PU")
highlight_PR <- TV_join %>% filter(treatment == "PR")

ggplot(data = TV_join, aes(x = Eri_vag, y = NDVI_avg)) +
  geom_point(size = 2)+
  geom_smooth(method = "lm", se = TRUE, color = "black")+
   geom_point(data=highlight_CT, aes(x=Eri_vag,y=NDVI_avg), size = 2, color = "#a6611a")+
  geom_point(data=highlight_EX, aes(x=Eri_vag,y=NDVI_avg), size = 2, color = "#dfc27d")+
  geom_point(data=highlight_PU, aes(x=Eri_vag,y=NDVI_avg), size = 2, color = "#80cdc1")+
  geom_point(data=highlight_PR, aes(x=Eri_vag,y=NDVI_avg), size = 2, color = "#018571")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Eri vag")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.3)+ #adds R2 
  #stat_cor(method = "pearson", label.x = .4)+ #adds p-value 
  #theme(legend.position="bottom")+
  facet_grid(site ~ .)

  #facet_grid(. ~ site)

ggsave("../Figures/Teamvole_regrestion_Erivag_NDVI.jpeg")

Rub cha vs NDVI

highlight_CT <- TV_join %>% filter(treatment == "CT")
highlight_EX <- TV_join %>% filter(treatment == "EX")
highlight_PU <- TV_join %>% filter(treatment == "PU")
highlight_PR <- TV_join %>% filter(treatment == "PR")

ggplot(data = TV_join, aes(x = Rub_cha, y = NDVI_avg)) +
  geom_point(size = 2)+
  geom_smooth(method = "lm", se = TRUE, color = "black")+
   geom_point(data=highlight_CT, aes(x=Rub_cha,y=NDVI_avg), size = 2, color = "#a6611a")+
  geom_point(data=highlight_EX, aes(x=Rub_cha,y=NDVI_avg), size = 2, color = "#dfc27d")+
  geom_point(data=highlight_PU, aes(x=Rub_cha,y=NDVI_avg), size = 2, color = "#80cdc1")+
  geom_point(data=highlight_PR, aes(x=Rub_cha,y=NDVI_avg), size = 2, color = "#018571")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Rub cha")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  #stat_regline_equation(label.x=0.2)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.2)+ #adds R2 
  #stat_cor(method = "pearson", label.x = .2)+ #adds p-value 
  #theme(legend.position="bottom")+
  facet_grid(site ~ .)

  #facet_grid(. ~ site)

ggsave("../Figures/Teamvole_regrestion_Rubcha_NDVI.jpeg")

LTER NDVI/COVER

Read in LTER NDVI data

LTER_NDVI_orig <- read.csv("../Data/2011_2018_LMAT_NDVI.csv", header = TRUE)

Prep the NDVI data

Subset by year, correct DOY column to reflect year

LTER_NDVI_11 <- LTER_NDVI_orig %>% filter(year == "2011") %>%
  mutate(doy_adj = doy)

LTER_NDVI_14 <- LTER_NDVI_orig %>% filter(year == "2014")%>%
  mutate(doy_adj = doy+(365*3))

LTER_NDVI_15 <- LTER_NDVI_orig %>% filter(year == "2015")%>%
  mutate(doy_adj = doy+(365*4))

LTER_NDVI_16 <- LTER_NDVI_orig %>% filter(year == "2016")%>%
  mutate(doy_adj = doy+(365*5))

LTER_NDVI_17 <- LTER_NDVI_orig %>% filter(year == "2017")%>%
  mutate(doy_adj = doy+(365*6))

LTER_NDVI_18 <- LTER_NDVI_orig %>% filter(year == "2018")%>%
  mutate(doy_adj = doy+(365*7))

merge data back together

LTER_NDVI_CLEAN           <- Reduce(MyMerge, list(LTER_NDVI_11, LTER_NDVI_14, LTER_NDVI_15, LTER_NDVI_16, LTER_NDVI_17, LTER_NDVI_18)) 

LTER_NDVI_CLEAN$time <- LTER_NDVI_CLEAN$time %>%
  recode("initial " = "initial")

Line graph of NDVI

#specify factor levels for fertilization and clipping treatments 
LTER_NDVI_CLEAN$plot <- factor(LTER_NDVI_CLEAN$plot, levels = c("F2", "F5", "F10", "CT"))


ggplot(data = LTER_NDVI_CLEAN, aes(x = doy_adj, y = avg_ndvi, color = plot)) +
  scale_color_manual(values = c("CT" = "#d8b365", 
                                "F2" = "#41b6c4",
                                "F5" = "#2c7fb8",
                                "F10" = "#253494"))+
  
  theme_light()+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1),
        axis.text.y = element_text(color = "black", size = 8),
        axis.line = element_line(colour = "black", linetype = "solid"), 
        legend.position="none",
        plot.title = element_text(hjust = .5))+
  geom_line(size = .75) +
  geom_point(size = 1.25) + 
  labs(y = "NDVI", x = "DOY")+ 
  theme(axis.title.y = element_text(face = "bold"))+
  theme(axis.title.x = element_text(face = "bold"))+
  scale_x_break(c(218, 1253)) +
  scale_x_break(c(1321, 1618)) +
  scale_x_break(c(1685, 1971)) +
  scale_x_break(c(2040, 2348)) +
  scale_x_break(c(2425, 2727))

ggsave("../Figures//LTER_NDVI_linegraph.jpeg")

## Saving 7 x 5 in image

Read in LTER cover data

LTER_cover_orig <- read.csv("../Derived CSVs/LTER_06MAT_relcov_funcgrp_allyrs.csv", header = TRUE)

Prep the cover data

Subset by year based on # of quads measured per plot in a given year

In 2008 we measured 4 quadrants/plot
In 2015 & 2020 we measured with 5 quadrants/plot
In all other years we measured with 8 quadrants

Q4 <- subset(LTER_cover_orig, year== "2008")
Q5 <- subset(LTER_cover_orig, year== "2015" | year== "2020")
Q8 <- subset(LTER_cover_orig, year== "2010" | year== "2011" | year== "2012" | year== "2013" | year== "2014" | year== "2016" | year== "2017" | year== "2018" | year== "2019" |year== "2010" |year== "2021")

Sum across quadrats in a plot

Q4_sum_avg_quad <- (Q4) %>% 
  group_by(year, block, plot, species) %>%
  summarise_at(vars(rel_cov), list(sum_quad = sum), na.rm = TRUE)

Q5_sum_avg_quad <- (Q5) %>% 
  group_by(year, block, plot, species)%>% 
summarise_at(vars(rel_cov), list(sum_quad = sum), na.rm = TRUE)

Q8_sum_avg_quad <- (Q8) %>% 
  group_by(year, block, plot, species) %>% 
  summarise_at(vars(rel_cov), list(sum_quad = sum), na.rm = TRUE)

Create functions/vectors for calculating average rel_cov bases on # of quadrats then fill a new column with the average

avg_Q4_quad <- c(Q4_sum_avg_quad$sum_quad/4)
avg_Q5_quad <- c(Q5_sum_avg_quad$sum_quad/5)
avg_Q8_quad <- c(Q8_sum_avg_quad$sum_quad/8)

#creates and fills a new column with the average (by diving that sum by the number of reps specified)
Q4_sum_avg_quad["avg_quad"] <-avg_Q4_quad
Q5_sum_avg_quad["avg_quad"] <-avg_Q5_quad
Q8_sum_avg_quad["avg_quad"] <-avg_Q8_quad

Merge all the years back together

MyMerge       <- function(x, y){
  df            <- merge(x, y, all = TRUE)
  rownames(df)  <- df$Row.names
  df$Row.names  <- NULL
  return(df)
}
Avg_quad           <- Reduce(MyMerge, list(Q4_sum_avg_quad, Q5_sum_avg_quad, Q8_sum_avg_quad))

Average across blocks

Avg_block <- (Avg_quad) %>% 
  group_by(year, plot, species) %>% 
  summarise_at(vars(avg_quad), list(avg_block = mean, block_sd = sd), na.rm = TRUE)

LTER_cover_CLEAN <- Avg_block %>%
  filter(species == "Bet nan" | species == "Eri vag" | species == "Rub cha") %>%
  select(-block_sd)

Join NDVI and COVER data

Join the dataset, drop columns we don’t need

NDVI_initial <- LTER_NDVI_CLEAN %>%
  filter(time == "initial")

LTER_initial <- left_join(NDVI_initial,LTER_cover_CLEAN, by= c("year", "plot")) %>%
  select(-doy_adj, -doy) %>%
  pivot_wider(names_from = species, values_from = avg_block) %>%
  rename("Bet_nan" = "Bet nan", "Eri_vag" = "Eri vag", "Rub_cha" = "Rub cha")

#######

NDVI_half <- LTER_NDVI_CLEAN %>%
  filter(time == "half")

LTER_half <- left_join(NDVI_half,LTER_cover_CLEAN, by= c("year", "plot")) %>%
  select(-doy_adj, -doy)%>%
  pivot_wider(names_from = species, values_from = avg_block) %>%
  rename("Bet_nan" = "Bet nan", "Eri_vag" = "Eri vag", "Rub_cha" = "Rub cha")

#######

NDVI_peak <- LTER_NDVI_CLEAN %>%
  filter(time == "peak")

LTER_peak <- left_join(NDVI_peak,LTER_cover_CLEAN, by= c("year", "plot")) %>%
  select(-doy_adj, -doy)%>%
  pivot_wider(names_from = species, values_from = avg_block) %>%
  rename("Bet_nan" = "Bet nan", "Eri_vag" = "Eri vag", "Rub_cha" = "Rub cha")

#######

NDVI_post <- LTER_NDVI_CLEAN %>%
  filter(time == "post-peak")

LTER_post <- left_join(NDVI_post,LTER_cover_CLEAN, by= c("year", "plot")) %>%
  select(-doy_adj, -doy)%>%
  pivot_wider(names_from = species, values_from = avg_block) %>%
  rename("Bet_nan" = "Bet nan", "Eri_vag" = "Eri vag", "Rub_cha" = "Rub cha")

Scatter plot of NDVI vs species cover

Bet nan vs NDVI

INITIAL

#specify factor levels for treatments
LTER_initial$plot <- factor(LTER_initial$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_in <- LTER_initial %>% filter(plot == "CT")
highlight_F2_in <- LTER_initial %>% filter(plot == "F2")
highlight_F5_in <- LTER_initial %>% filter(plot == "F5")
highlight_F10_in <- LTER_initial %>% filter(plot == "F10")

Bet_nan_initial <-  ggplot(data = LTER_initial, aes(x = Bet_nan, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_in, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_in, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_in, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_in, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Bet nan")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  #stat_cor(method = "pearson", label.x = .3)+ #adds p-value 
  theme(legend.position="none")


print(Bet_nan_initial)

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

HALF

#specify factor levels for treatments
LTER_half$plot <- factor(LTER_half$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_half <- LTER_half %>% filter(plot == "CT")
highlight_F2_half <- LTER_half %>% filter(plot == "F2")
highlight_F5_half <- LTER_half %>% filter(plot == "F5")
highlight_F10_half <- LTER_half %>% filter(plot == "F10")

Bet_nan_half <-  ggplot(data = LTER_half, aes(x = Bet_nan, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_half, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_half, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_half, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_half, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Bet nan")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  #stat_cor(method = "pearson", label.x = .3)+ #adds p-value 
  theme(legend.position="none")


print(Bet_nan_half)

PEAK

#specify factor levels for treatments
LTER_peak$plot <- factor(LTER_peak$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_peak <- LTER_peak %>% filter(plot == "CT")
highlight_F2_peak <- LTER_peak %>% filter(plot == "F2")
highlight_F5_peak <- LTER_peak %>% filter(plot == "F5")
highlight_F10_peak <- LTER_peak %>% filter(plot == "F10")


Bet_nan_peak <-  ggplot(data = LTER_peak, aes(x = Bet_nan, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_peak, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_peak, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_peak, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_peak, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Bet nan")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  #stat_cor(method = "pearson", label.x = .3)+ #adds p-value 
  theme(legend.position="none")


print(Bet_nan_peak)

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

POST-PEAK

#specify factor levels for treatments
LTER_post$plot <- factor(LTER_post$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_post <- LTER_post %>% filter(plot == "CT")
highlight_F2_post <- LTER_post %>% filter(plot == "F2")
highlight_F5_post <- LTER_post %>% filter(plot == "F5")
highlight_F10_post <- LTER_post %>% filter(plot == "F10")


Bet_nan_post <-  ggplot(data = LTER_post, aes(x = Bet_nan, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_post, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_post, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_post, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_post, aes(x=Bet_nan,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Bet nan")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  #stat_cor(method = "pearson", label.x = .3)+ #adds p-value 
  theme(legend.position="none")


print(Bet_nan_post)

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

Combine the graphs

(Bet_nan_initial+Bet_nan_half)/(Bet_nan_peak+Bet_nan_post) + plot_annotation(tag_levels = 'A')

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

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

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

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

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

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

ggsave("../Figures/LTER_Betnan_NDVI_reg.jpg")

## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

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

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

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

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

Eri vag vs NDVI

INITIAL

#specify factor levels for treatments
LTER_initial$plot <- factor(LTER_initial$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_in <- LTER_initial %>% filter(plot == "CT")
highlight_F2_in <- LTER_initial %>% filter(plot == "F2")
highlight_F5_in <- LTER_initial %>% filter(plot == "F5")
highlight_F10_in <- LTER_initial %>% filter(plot == "F10")


Eri_vag_in <-  ggplot(data = LTER_initial, aes(x = Eri_vag, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_in, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_in, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_in, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_in, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Eri vag")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  stat_cor(method = "pearson", label.x = .15)+ #adds p-value 
  theme(legend.position="none")


print(Eri_vag_in)

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing non-finite values (stat_cor).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

HALF

#specify factor levels for treatments
LTER_half$plot <- factor(LTER_half$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_half <- LTER_half %>% filter(plot == "CT")
highlight_F2_half <- LTER_half %>% filter(plot == "F2")
highlight_F5_half <- LTER_half %>% filter(plot == "F5")
highlight_F10_half <- LTER_half %>% filter(plot == "F10")


Eri_vag_half <-  ggplot(data = LTER_half, aes(x = Eri_vag, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_half, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_half, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_half, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_half, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Eri vag")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  stat_cor(method = "pearson", label.x = .15)+ #adds p-value 
  theme(legend.position="none")


print(Eri_vag_half)

PEAK

#specify factor levels for treatments
LTER_peak$plot <- factor(LTER_peak$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_peak <- LTER_peak %>% filter(plot == "CT")
highlight_F2_peak <- LTER_peak %>% filter(plot == "F2")
highlight_F5_peak <- LTER_peak %>% filter(plot == "F5")
highlight_F10_peak <- LTER_peak %>% filter(plot == "F10")


Eri_vag_peak <-  ggplot(data = LTER_peak, aes(x = Eri_vag, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_peak, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_peak, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_peak, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_peak, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Eri vag")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  stat_cor(method = "pearson", label.x = .15)+ #adds p-value 
  theme(legend.position="none")


print(Eri_vag_peak)

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing non-finite values (stat_cor).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

POST-PEAK

#specify factor levels for treatments
LTER_post$plot <- factor(LTER_post$plot, levels = c("CT","F2","F5","F10"))

highlight_CT_post <- LTER_post %>% filter(plot == "CT")
highlight_F2_post <- LTER_post %>% filter(plot == "F2")
highlight_F5_post <- LTER_post %>% filter(plot == "F5")
highlight_F10_post <- LTER_post %>% filter(plot == "F10")


Eri_vag_post <-  ggplot(data = LTER_post, aes(x = Eri_vag, y = avg_ndvi)) +
  geom_point(size = 2)+
    geom_smooth(method = "lm", se = TRUE, color = "black")+
  geom_point(data=highlight_CT_post, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#d8b365")+
  geom_point(data=highlight_F2_post, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#41b6c4")+
  geom_point(data=highlight_F5_post, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#2c7fb8")+
  geom_point(data=highlight_F10_post, aes(x=Eri_vag,y=avg_ndvi), size = 2, color = "#253494")+
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(color = "black", size = 8, angle = 45, hjust = 1))+
  theme(axis.text.y = element_text(color = "black", size = 8))+
  theme( axis.line = element_line(colour = "black", 
                      size = 1, linetype = "solid"))+
  labs(y = "NDVI", x = "Eri vag")+ 
  theme(axis.title.y = element_text(face = "bold"))+
    theme(axis.title.x = element_text(face = "bold"))+
  ylim(0.2, 0.9)+
  #stat_regline_equation(label.x=0.3)+ #adds regression line equation
  #stat_cor(aes(label=..rr.label..), label.x=0.35)+ #adds R2 
  stat_cor(method = "pearson", label.x = .15)+ #adds p-value 
  theme(legend.position="none")


print(Eri_vag_post)

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing non-finite values (stat_cor).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

Combine the graphs

(Eri_vag_in+Eri_vag_half)/(Eri_vag_peak+Eri_vag_post) + plot_annotation(tag_levels = 'A')

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

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing non-finite values (stat_cor).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

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

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing non-finite values (stat_cor).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

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

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing non-finite values (stat_cor).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

ggsave("../Figures/LTER_Erivag_NDVI_reg.jpg")

## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing non-finite values (stat_cor).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

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

## Warning: Removed 2 rows containing non-finite values (stat_smooth).

## Warning: Removed 2 rows containing non-finite values (stat_cor).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

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

## Warning: Removed 1 rows containing non-finite values (stat_smooth).

## Warning: Removed 1 rows containing non-finite values (stat_cor).

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

NDVI_exploration

Nicole

3/21/2022

Team vole NDVI/COVER

Load in NDVI data

Prep the Team Vole NDVI data

Box plot - NDVI

Load in Team Vole cover data

Prep TV cover

Join NDVI and COVER data

Scatter plot of NDVI vs species cover

LTER NDVI/COVER

Read in LTER NDVI data

Prep the NDVI data

Line graph of NDVI

Read in LTER cover data

Prep the cover data

Join NDVI and COVER data

Scatter plot of NDVI vs species cover

Bet nan vs NDVI

Eri vag vs NDVI